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A  GEOMETRIC  OPTICS  MODEL  FOR  CALCULATING  THS  FELD  STRENGTH  OF  ELECTROMAG¬ 
NETIC  WAVES  IN  THE  PRESENCE  OF  A  TROPOSPHERIC  DUCT 

Raymond  ?.  Hasky 
University  of  Dayton,  1977 

Major  Professor*  Dr.  3.  K.  Schmidt 

The  theory  and  development  of  a  geometric  optics  model  for  analyzing 
the  effects  of  anomalous  atmospheric  refraction  on  the  field  strength  of 
radio  frequency  emitters  is  presented.  The  model  is  derived  from  Fermat's 
principle  which  defines  the  Euler-Lagrange  equations  of  a  ray  and  the  ray 
optical  path  length.  This  model  is  applicable  to  radio  propagation  above 
30  MHz  where  ionospheric  effects  are  generally  negligible.  Given  an  iso¬ 
tropic  emitter  of  known  frequency,  polarization,  pulse  width,  and  alti¬ 
tude,  the  free  space  normalized  power  density  and  field  strength  are  cal¬ 
culated  as  a  function  of  altitude  and  distance  along  the  earth’s  surface. 
Computations  are  made  by  constructing  a  ray  trace  and  utilizing  the  Jacob¬ 
ian  of  the  transformation  between  ray  space  coordinates  and  wavefront 
surface  coordinates  to  solve  for  the  time-averaged  free  space  normalized 
roynting  vector. 

The  atmosphere  is  treated  as  a  non-raagnetic  Inhomogeneous  medium 
which  is  othenrise  linear  (|i  and  €  are  independent  of  the  fields)  and 
isotropic  ((.!  and  5  are  scalars).  The  atmospheric  refractive  index  is 
modeled  as  a  function  of  altitude  above  a  spherical  earth  surface  by 
either  a  standard  Central  Radio  Propagation  Laboratory  (CRFL)  exponential 
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function  or  by  measured  profiles  obtained  from  meteorological  soundings 
(horizontal  homogeneity  in  the  refractive  index  may  be  assumed  over 
relatively  large  distances).  Fields  that  are  incident  upon  the  earth's 
surface  are  specularly  reflected  and  attenuated  by  a  Fresnel  reflection 
coefficient  and  a  surface  roughness  factor  which  is  dependent  upon  the 
angle  of  incidence  and  standard  deviation  of  the  surface  about  a  mean 
height. 

Results  obtained  with  this  model  are  compared  to  experimental 
waveguide  mode  theory  field  strength  calculations  for  a  low  altitude 
superrefractive  atmospheric  layer,  or  tropospheric  duct,  lying  along  a 
230  nautical  mile  path  off  the  California  coast  between  San  Diego  and 
Guadalupe  Island.  Field  strength  calculations  are  presented  at  65,  17C, 
520,  and  3300  KHz  first  using  the  average  measured  refractive  index  pro¬ 
file  for  the  Guadalupe  Island  duct  and  then  a  trilinear  approximation  to 
that  profile.  The  results  show  a  fair  amount  of  agreement  with  the  ex¬ 
perimental  and  waveguide  data  at  the  lower  frequencies,  with  good  to 
excellent  comparisons  at  520  and  3300  KHz  especially  when  the  Guadalupe 
Island  refractive  index  profile  is  used.  Calculations  are  restricted 
to  the  region  within  the  duct  because  of  limitations  in  the  geometric 
optics  method,  which  are  discussed  at  length. 
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CHAPTER  I 


INTRODUCTION  AND  DISCUSSION  OF  THE  THESIS 

Background 

The  beginning  of  large  scale  research  in  electromagnetic  wave  prop¬ 
agation  over  the  surface  of  the  earth  was  prompted  by  Marconi's  demon¬ 
stration  in  1901  that  signals  could  be  transmitted  across  the  Atlantic 
Ocean.  Much  of  the  work  concentrated  on  providing  a  description  of  four 
mechanisms  that  had  been  proposed  to  account  for  over-the-horizon  trans¬ 
mission!  diffraction  by  the  earth,  reflection  from  am  elevated  layer  of 
ionized  gases  (i.e.,  the  ionosphere),  atmospheric  refraction,  and  sur¬ 
face  waves  at  the  boundary  between  differing  dielectric  media.1  Papers 
by  Zenneck,  Sommerfeld,  and  Watson  concluded  that  long  distance  propaga¬ 
tion  up  to  30  MHz  was  explainable  in  terms  of  all  of  these  mechanisms 
except  atmospheric  refraction,  which  was  found  to  have  a  negligible  ef¬ 
fect.1*  2 

During  the  period  immediately  preceding  and  following  World  War  II, 
the  useable  radio  spectrum  was  extended  from  30  MHz  to  approximately 
300  GHz.  In  this  region  the  effects  of  surface  waves,  earth  diffraction, 
and  ionospheric  reflection  are  generally  absent,  while  refraction  by  the 
lower  atmosphere  or  troposphere,  becomes  the  central  mechanism  which  en- 
ables  long  distance  signal  transmission.  As  will  be  shown  in  Chapters 
II  and  IV,  it  is  not  the  absolute  value  of  the  refractive  index  of  air 
that  is  important  in  describing  radio  propagation  at  these  wavelengths, 
but  rather  its  rate  of  change. 


1 


Under  "standard"  conditions  the  atmospheric  index  of  refraction 

remains  essentially  constant  in  horizontal  directions  while  decreasing 

1  4 

exponentially  with  increasing  altitude.  ’  Srom  Snell’s  law  of  op¬ 
tics  (Chapter  III)  this  negative  gradient  in  the  refractive  index  nor¬ 
mally  sends  horizontally  launched  radio  waves  around  the  curvature  of 
the  earth  so  that  the  radio  horizon  is  about  4/3  the  distance  to  the 
geometric  horizon.  If  the  gradient  has  a  greater  than  normal  magnitude 
the  waves  will  be  refracted  even  further  around  the  earth's  surface, 
thus  extending  the  radar  horizon,  '..‘hen  the  gradient  becomes  suffici¬ 
ently  large,  the  waves  may  have  the  same  curvature  as  that  of  the  earth 
and  follow  the  earth's  surface  indefinitely,  a  condition  referred  to  as 
superrefraction,  trapping,  or  ducting.  Conversely,  a  smaller  than  usual 
change  in  the  refractive  index  results  in  substandard  propagation  or 
subrefraction. ^ 

The  refractive  index  of  air  at  radio  frequencies  is  a  function  of 
total  atmospheric  pressure,  water  vapor  content,  and  temperature  (Chap¬ 
ter  II)  and  is  therefore  subject  to  local  meteorological  conditions. 
Greater  than  normal  refraction  generally  occurs  when  the  temperature  in¬ 
creases  or  the  humidity  decreases  rapidly  with  height,  such  as  when  a 
warm  dry  air  mass  from  the  land  is  blown  out  over  a  cooler  sea  surface, 
or  when  heat  is  radiated  from  the  earth's  surface  at  night  and  the  ground 
is  moist.  Such  conditions  are  called  temperature  and  humidity  inver¬ 
sions,  with  humidity  inversions  usually  being  the  dominant  factor  in  de¬ 
termining  atmospheric  refraction  properties.  As  a  result,  superrefract- 

ive  and  nearly  superrefractive  propagating  conditions  are  more  prevalent 
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over  oceans  than  over  land  surfaces,  ’  ’  Superrefractive  conditions, 

’••'ben  they  do  exist  over  land,  tend  to  exhibit  more  pronounced  daily  var¬ 
iations  than  over  the  sea,  since  land  masses  change  temperature  more 


3 

rapidly  than  do  sea  bodies.  In  either  case,  superrefraction  occurs  most 
frequently  under  fair  weather  conditions  when  the  atmosphere  is  verti¬ 
cally  stratified  and  there  is  little  or  no  air  mixing  due  to  turbulence. 
Propagation  is  more  nearly  normal  when  the  atmosphere  is  cold  or  stormy 
or  when  very  rough  terrain  and  high  winds  contribute  to  the  mixing  proc¬ 
ess.^ 

Long  distance  propagation  due  to  superrefraction  is  possible  with 

only  a  relatively  small  change  in  the  refractive  index.  A  mere  gradient 

of  (-1.57)l0~ ^  parts  per  meter  in  altitude  for  the  index  will  result  in 

7 

ducting.  However,  as  Skolnik  illustrates,  the  results  may  be  dramatic, 

as  in  the  case  of  a  ground  based  200  KHz  radar  in  Bombay,  India,  which 

frequently  received  target  echoes  from  points  in  Arabia  at  ranges  of 

1000  to  1500  miles  during  World  War  II.  Such  extremely  long  distance 

operation  is  possible  since  the  atmospheric  duct  behaves  like  a  leaky 

waveguide,  with  the  ground  plus  the  thermal  or  humidity  inversion  layers 

acting  as  lossy  waveguide  walls.  While  most  ducts  are  bounded  by  the 

earth  on  the  bottom  and  an  upper  inversion  layer  usually  not  more  than 

several  hundred  feet  high,  elevated  earth-detached  ducts  have  been  re- 
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ported  up  to  altitudes  of  several  thousand  feet.' *  ’ 

7 

As  an  example  of  subrefractive  conditions,  Skolnik  cites  an  in¬ 
stance  when  radars  off  Fisher's  Island,  New  York,  were  unable  to  see 
31ock  Island  only  22  miles  away,  although  it  was  optically  visible.  This 
occurred  because  the  refractive  index  gradient,  rather  than  bending  the 
signals  around  the  earth's  curvature,  refracted  the  signals  away  from 
the  earth,  thus  reducing  the  range  to  the  radar  horizon.  Under  certain 
circumstances  subrefraction  may  result  from  a  low  flying  fog  for  which 
part  of  the  water  in  the  air  changes  from  a  gaseous  to  a  liquid  state. 
Although  the  total  amount  of  water  in  the  air  remains  unchanged,  only  the 


water  in  a  vapor  form  contributes  significantly  to  the  air's  refractive 

index,  making  the  index  lower  than  normal  at  the  surface,  and  conse- 

7 

quently  yielding  substandard  propagation  conditions.' 

Statement  of  the  Problem 

i'hile  surveying  the  literature  for  the  effects  of  refraction  on 
propagation  at  radio  frequencies,  it  becomes  evident  that  the  vast  ma¬ 
jority  of  analytical  work  in  this  area  has  been  conducted  by  using  prop¬ 
agation  models  that  are  either  extremely  simplistic  or  exceedingly  com¬ 
plex  and  specialized.  The  simpler  approaches  generally  assume  some  form 
of  ’'standard"  vertical  refractive  index  profile^*  ^  such  as  a  constant 
gradient  index  or  "4/3  effective  earth  radius"  model  (Chapter  II),  which 
is  not  representative  of  superrefractive  or  subrefractive  atmospheric 

conditions.  The  more  complex  methods,  which  use  waveguide  mode  theory 
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to  solve  for  tropospheric  refraction,  *  J  are  generally  limited  to 
atmospheric  structures  having  only  one  or  two  inversion  layers  since  the 
difficulty  of  solution  increases  rapidly  with  the  number  of  changes  in 
the  vertical  profile.  In  either  case,  measured  refractive  index  data 
is  usually  replaced  by  approximate  profile  models  which,  as  Chapter  VI 
indicates,  may  have  a  noticeable  effect  on  results. 

The  purpose  of  this  thesis  is  to  present  the  theory  and  development 
of  an  accurate  yet  practical  model  for  analyzing  the  effects  of  atmos¬ 
pheric  refraction  on  the  field  strength  of  radio  frequency  emitters  in 
the  30  to  100  GHz  spectrum.  Results  obtained  from  this  model  are 
compared  with  theoretical  and  experimental  data  recently  reported  by  Pap- 

14 

pert  and  Goodhart  for  a  ground-based  tropospheric  duct  lying  along  a 
280  nautical  mile  path  between  San  Diego  and  Guadalupe  Island.  Field 
strength  calculations  are  made  at  65,  170,  520,  and  3300  ?3iz  using  the 
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measured  refractive  index  profile  and  a  trilinear  approximation  to  that 
profile. 

Since  this  is  a  study  of  refraction  effects,  no  attempt  is  made  to 
include  signal  attenuation  due  to  precipitation  scattering  or  gaseous 
absorption,  both  of  which  are  widely  discussed  in  the  literature.  Fur¬ 
thermore,  the  choice  of  frequency  range  in  this  study  eliminates  the 
problem  of  partial  penetration  and  reflection  from  ionized  layers  in  the 
atmosphere.  Consequently,  the  atmosphere  is  treated  as  an  isotropic  and 
linear  medium  (i.e.,  the  permeability  and  permittivity  are  scalar  quan¬ 
tities  that  are  independent  of  the  fields  in  the  medium)  which  would  not 
be  the  case  if  ionospheric  effects  were  present. 

Method  of  Approach 

One  of  the  major  concerns  that  arose  early  in  this  study  regarded 
the  selection  of  a, suitable  approach  for  modeling  propagation  through  a 
refractive  medium.  Classically  the  solutions  for  radio  propagation 
problems  fall  into  two  general  categories i  geometric  optics  or  ray  the¬ 
ory,  and  physical  optics  or  wave  theory.  Each  has  its  advantages  and 
limitations,  depending  on  the  type  of  problem  being  addressed. 

Geometric  optics  is  usually  a  simpler  approach  than  physical  optics 
since  it  describes  the  propagation  of  waves  along  rays  according  to  ele¬ 
mentary  geometric  laws  without  regard  to  wavelengths  or  phases.  The  rays 
are  defined  as  normals  to  the  surfaces  of  constant  phase  of  the  wave- 
front.  ”hen  the  medium  is  isotropic,  the  index  of  refraction  may  be  re¬ 
garded  to  be  a  real  quantity,  and  the  rays  are  found  to  lie  along  the 
direction  of  energy  propagation.  For  anisotropic  media,  however,  the  re¬ 
fractive  index  becomes  complex  and  the  rays  may  not  follow  the  direction 
of  energy  flow.1^ 
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The  attractiveness  of  geometric  optics  is  that  it  retains  its  sim¬ 
plicity  even  for  complicated  refractive  index  structures,  while  genera¬ 
ting  easily  interpreted  ray  patterns  which  show  the  effect  of  refraction 
on  an  emitter' s  radiation  pattern.  However  its  main  drawbacks  are  that 
it  is  not  valid  for  diffraction  or  complex  scattering  problems,  nor  is 
it  suitable  in  regions  of  rapid  refractive  index  changes  (within  a  wave¬ 
length  of  distance)  or  rapid  ray  divergence,  such  as  near  sources  or 
focal  points.  Any  use  of  ray  optics  must  therefore  include  careful  con¬ 
sideration  of  its  limitations,  with  the  knowledge  that  it  is  only  a  lim¬ 
iting  form  of  physical  optics.1'  1^ 

Physical  optics,  by  contrast,  results  from  a  solution  of  the  wave 
equation,  and  introduces  wavelengths,  phases,  and  mode  concepts.  It  is 
valid  for  all  diffraction  and  scattering  problems  and  is  unaffected  by 
rapid  changes  in  the  refractive  index  or  by  the  presence  of  sources  or 
focal  points.1^  Unfortunately  this  method  rapidly  becomes  unmanageable 
for  complicated  atmospheric  structures  to  the  point  where  it  is  nearly 
useless.  Furthermore,  while  presenting  "exact"  solutions  to  propagation 
problems,  its  accuracy  generally  exceeds  that  of  available  measured  re- 
fractivity  data.  Even  for  complex  scattering  problems  it  becomes  so 

difficult  an  approach  that  many  analytic  methods  still  rely  on  the  ap- 

1  17 

proximations  of  ray  theory  whenever  feasible.  ' 

Considering  the  applicability  of  ray  versus  wave  theory,  the  former 
was  selected  for  use  in  this  study  because  of  its  inherent  simplicity. 
Since  earth  diffraction  is  not  regarded  as  a  significant  contributor  to 
long  range  propagation  at  the  wavelengths  being  considered,  it  may  be 
eliminated  as  an  obstacle  to  the  use  of  geometric  optics.  Scattering  in 
this  paper  is  limited  to  specular  reflection  from  the  earth,  since  dif¬ 
fuse  reflection  becomes  an  extremely  complicated  function  of  incidence 
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.ir.^le,  wavelength,  polarization,  surface  roughness,  and  surface  electri¬ 
cal  properties.  The  problem  of  modeling  diffuse  scattering  from  rough 
surfaces  is  so  difficult  that  many  authors  choose  to  describe  it  in 
terms  of  a  summation  of  specular  reflections  from  randomly  oriented 
olane  facet  surfaces. Specular  reflection  on  the  other  hand,  may 
be  included  in  the  geometric  optics  approach  by  assigning  complex  spec¬ 
ular  scattering  coefficients  to  the  reflected  signal,  as  shown  in  Chap¬ 
ter  V. 

numerous  geometric  optics  models  exist  which  use  a  variety  of  tech¬ 
niques  to  compute  field  strength  from  the  ray  trajectories  of  a  refract- 
i  q_22 

ed  field.  Each  of  these,  however,  has  considerable  limitations 

either  in  the  selection  of  atmospheric  models,  method  of  field  strength 
calculation,  absence  of  a  reasonable  earth  reflection  model,  or  general 
lack  of  versatility.  Consequently  it  was  necessary  to  develop  a  new  com¬ 
puter  program  for  this  thesis  (Chapter  V)  which  solves  the  ray  trajectory 
equations  derived  in  Chapter  III  and  calculates  emitter  field  strength  as 
a  function  of  height  and  distance  along  the  earth's  surface.  While  this 
simulation  represents  original  work,  portions  of  the  field  strength  and 
numerical  integration  algorithms  were  obtained  from  a  computer  program 
described  in  Reference  21. 


CHAPTER  II 


RADIO  REFRACTIVTTY  OF  AIR 


As  previously  mentioned,  long  range  propagation  at  radio  frequen¬ 
cies  above  30  KHz  is  mainly  attributable  to  vertical  changes  in  the  re¬ 
fractive  index  of  air.  Refractive  index  is  defined  as  the  ratio  of  the 
velocity  of  light  in  free  space  c  to  the  phase  velocity  v^  of  a  field 
traveling  through  a  given  medium.  Thus  the  refractive  index  of  air  is 


.y/ji* 

Vv  o 


where  |i  and  €  are  the  permeability  and  permittivity  of  air,  and  |iQ,  €q 
are  the  respective  quantities  in  free  space. 

For  good  dielectric  and  non-magnetic  media  such  as  air 

4  -  404r  *  40  (2) 

€  «  €  €  (3) 

o  r 

where  |ir  and  €r  are  the  relative  permeability  and  permittivity.  Eq.  (l) 
then  may  be  written 
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Since  n  w  1  for  air,  Eq.  (4)  may  be  approximated  by  the  binomial  series 


, ,  ,  x  x  , 
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such  that 
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n  =  V  1  +  (€r  -  l)  *  1  +  - 

3a.  (6)  is  generally  rewritten  to  define  the  refractivity  «  of  air  as 


(6) 


l<  =  (n  -  1)-10°  (7) 

3he  use  of  3  is  widespread  in  radio  meteorology  since  it  reflects  the 
difference  between  the  refractive  indices  of  air  and  vacuum  in  units 
that  range  from  zero  to  several  hundreds.  At  sea  level  the  value  of  n 
is  approximately  n  =  1.00031,  which  corresponds  to  an  N  value  of  310. 

Other  commonly  used  units  include  the  "modified  indices"  given  by 


3  -  (n  -  1 


(8) 


and 

K  -  (n  -  1  +  (9) 

where  h  is  height  above  the  earth's  surface  and  a  is  the  radius  of  the 
earth.  As  will  be  shown  in  Chapter  IV,  the  radius  of  curvature  of  a 
nearly  horizontal  traveling  wave  is  inversely  proportional  to  the  verti¬ 
cal  gradient  dn/dh  of  the  air's  refractive  index.  Since  -l/4a  is  close 
to  the  observed  gradient  of  n  at  low  altitudes  under  standard  atmospheric 
conditions,^  the  3-unit  serves  to  eliminate  this  standard  decrease  of  N 
by  adding  the  quantity  (h/4a)*10^  to  if.  Onus  a  positive,  zero,  or  nega¬ 
tive  valued  3-gradient  represents  below  standard,  standard,  or  above 
standard  refractive  conditions,  respectively.  For  ducting,  when  a  sig¬ 
nal’s  radius  of  curvature  is  equal  to  that  of  the  earth,  dn/dh  =  -l/a. 

6 

Adding  (h/a)*10'  to  N  yields  the  K  unit,  so  that  M-gradients  are  always 
positive  except  for  ducting  conditions,  when  the  gradients  assume  zero  or 
negative  values.  The  major  advantage  of  using  refractivity,  or  IV-units, 
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lies  in  the  fact  that  N  is  not  modified  by  any  altitude  term  and  there¬ 
fore  reflects  the  true  distribution  of  the  refractive  index  n  with  re¬ 
spect  to  altitude.  The  gradient  of  N  is  generally  negative  even  for 
some  conditions  of  subrefraction.  For  ducting  the  gradient  must  sat¬ 
isfy 

gjj  *  -  j  -  -157  N-units/km  (10) 

The  refractivity  of  air  is  a  function  of  atmospheric  temperature, 
pressure,  and  water  vapor  content.  While  N  varies  with  frequency  in 
general,  it  is  essentially  frequency  independent  over  the  radio  spectrum 
being  considered  (30  KHz  to  100  GHz).  Numerous  studies  have  been  con¬ 
ducted  to  accurately  describe  the  radio  refractivity  of  air  in  terms  of 

atmospheric  properties.  The  generally  used  expression  for  N  comes  from 
'  24 

the  work  by  Smith  and  Weintraub  and  is  given  by 

p 

N  «  77.6  ^  +  72  |  +  (3.75)-105  ~  (11) 

where  is  the  pressure  of  dry  air  in  millibars  (mbar),  T  is  the  tem¬ 
perature  in  degrees  Kelvin  (°K),  and  e  is  the  partial  pressure  of  water 
vapor  in  mbar.  According  to  Smith  and  Weintraub,  Eq.  (ll)  has  an  over¬ 
all  accuracy  of  i  0.5  percent  of  N  for  frequencies  up  to  30  GHz  (exclud¬ 
ing  refractivity  dispersion  at  the  22  GHz  water  vapor  resonance  lines) 
and  for  normally  encountered  ranges  of  pressure,  temperature,  and  humid¬ 
ity.  This  accuracy  also  includes  errors  due  to  rounding  the  constants 
to  three  significant  figures.  The  constants  of  Eq.  (ll)  are  in  basic 
agreement  with  the  results  of  many  experimenters  whose  work  extends  up 
to  approximately  100  GHz.23 

Under  most  conditions  the  total  atmospheric  pressure  P  is 

P  «  P  +  p  +  e  «  F ,  +  e 
c  d  d 


(12) 


11 


Eq.  (11)  may  be 
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where  P  is  the  partial  pressure  of  carbon  dioxide, 
c 

written 

N  «  77.6  |  -  5.6  |  +  (3.75)-105  ~  (13) 

which  for  T  *  273  °K  («  0°  C)  may  be  simplified  to 

N  -  77.6  |  +  (3.73).105  ~  (14) 

or 

N  -  ZZi£  (P  +  4810  |)  (15) 

Ideally  the  refractive  index  is  directly  measurable  by  using  radio 
refractometers  which  are  sensitive  to  the  velocity  of  propagation  through 
the  air.  However  refractometers  are  expensive  and  often  difficult  to 
maintain  and  are  therefore  not  in  general  use.  The  bulk  of  refractive 
index  measurements  are  made  indirectly  from  meteorological  observations 
of  atmospheric  pressure,  temperature,  and  humidity,  which  are  then  con¬ 
verted  to  units  of  refractivity  through  Eqs.  (14)  or  (15).  These  meas¬ 
urements  are  typically  made  by  radiosondes  which  are  sent  aloft  on 
weather  balloons  at  thousands  of  weather  stations  around  the  world,  usu¬ 
ally  at  six-  or  twelve-hour  intervals. 

Radiosondes,  which  carry  sensors  that  detect  changes  in  pressure, 
temperature,  and  humidity,  transmit  their  data  to  the  ground  as  they 
rise  with  what  is  usually  assumed  to  be  a  constant  rate  of  ascent.  Un¬ 
fortunately  the  radiosonde  sensors  suffer  from  response  lags  and  some¬ 
times  rapid  component  aging,  which  can  adversely  affect  the  accuracy  of 
their  data.  Furthermore,  the  problem  of  refractivity  measurement  is 
complicated  by  the  fact  that  three  sets  of  data  are  required,  each  of 
which  contains  its  own  sensor  errors.  Thus  the  errors  in  the  data 
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generally  outweigh  any  discrepancies  in  the  constants  used  in  Eqs.  (14) 
or  (15). 23 

Numerous  models  of  "standard"  or  mean  atmospheres  exist.  They 
range  in  complexity  from  the  often  used  "4/3  earth  radius"  model  to  ex¬ 
ponential  and  bi-exponential  models.  Each  has  its  inherent  strengths 
and  weaknesses  and  its  range  of  applications.  For  example,  the  4/3 
earth  radius  model  (named  because  it  uses  an  effective  earth  whose  ra¬ 
dius  is  4/3  the  normal  earth  radius)  assumes  an  atmosphere  of  constant 
N-gradient  with  respect  to  altitude.  This  model  is  convenient  because 
of  its  computational  simplicity  and  because  it  allows  signals  to  propa¬ 
gate  along  straight  instead  of  curved  lines.  However  a  major  disadvan¬ 
tage  is  that  the  model  becomes  highly  inaccurate  at  altitudes  above  one 
or  two  kilometers,  a  fact  that  is  often  overlooked  by  many  of  its  users. 

The  model  which  most  generally  and  accurately  represents  standard 
atmospheric  refractivity  is  the  exponential  function 

N  -  N.  e  ~£ce(h  ’  *s>3  (16) 

where  Ng  is  the  refractivity  at  the  earth's  surface,  cg  is  the  exponen¬ 
tial  decay  rate  of  N,  h  is  height  above  sea  level,  and  h  is  the  surface 

s 

altitude  above  sea  level  at  which  N  is  measured.  This  model  has  been 

s 

adopted  for  use  by  the  National  Bureau  of  Standards  Central  Radio  Propa¬ 
gation  Laboratory  (CRPL)  along  with  a  table  of  values  for  N  and  c 

s  e 

This  model  and  its  associated  constants  from  the  CRH,  Exponential  Refer¬ 
ence  Atmosphere  (1958)  is  in  agreement  with  data  collected  by  the  Rocket 
?anel  and  the  Air  Force  Research  and  Development  Command.  Using  the 
mean  values  of  If,  c  ,  and  h  for  the  United  States,  Eq.  (16)  becomes 

H  -  3i3e"(°,l439h)  (17) 

where  h  is  in  kilometers.  Eq.  (17)  is  the  model  used  to  represent  a 


13 

standard  reference  atmosphere  in  this  paper.  A  variant  of  2q.  (17)  is 
the  bi-exponential  model  which  has  a  humidity  dependent  term  and  a  pres¬ 
sure  and  temperature  dependent  term.  The  bi-exponential  model  while 
providing  some  additional  modeling  flexibility  especially  in  regions  of 
significant  humidity  changes,  does  not  differ  appreciably  from  the 
single  term  model  for  most  general  cases.  J  In  any  event,  Eq.  (17) 
merely  serves  as  a  useful  analytic  standard  since  it  cannot  represent 
the  structure  of  N  for  any  given  location,  season,  or  time  of  day.  Thus 
to  realistically  study  the  effects  of  atmospheric  refraction  on  an  emit¬ 
ter's  field  distribution  for  a  given  set  of  meteorological  conditions, 

28 

it  becomes  necessary  to  use  measured  N  profiles. 

Finally,  the  discussion  of  refractivity  structure  has  centered  upon 
the  relationship  of  N  to  altitude,  and  has  neglected  changes  in  the  hor¬ 
izontal  direction.  An  examination  of  climatic  data  indicates  that  this 

is  reasonable  since  the  horizontal  changes  in  refractivity  are  slow  as 

23 

compared  to  the  vertical  changes.  Bean  indicates  that  it  would  be 
necessary  to  compare  sea  level  meteorological  stations  located  500  kil¬ 
ometers  apart  in  order  to  detect  a  difference  in  refractivity  values 
that  would  be  comparable  to  ascending  merely  one  kilometer  above  the 
location  of  either  station.  There  are,  as  Bean  reports,  several  special 
cases  where  the  horizontal  gradient  of  N  may  become  large,  and  these 
must  be  treated  on  an  individual  basis.  However  the  data  available  in¬ 
dicates  that  the  effect  of  horizontal  changes  in  N  is  generally  small 
and  can  usually  be  omitted  from  consideration. 


CHAPTER  III 


THEORY  OF  GEOMETRIC  OPTICS 

The  most  common  approach  to  developing  the  theory  of  geometric  op¬ 
tics  is  from  the  use  of  Snell's  law.  However,  considerable  insight  may 
be  gained  by  following  Kelso' s^  derivation  of  the  differential  ray 
equations  (Euler-Lagrange  equations)  from  Fermat’s  principle  in  optics. 
This  approach  directly  presents  solutions  for  the  optical  path  length 
of  a  ray,  the  time  of  propagation  along  a  ray,  and  the  relationship  of 
ray  height  to  its  projected  distance  along  the  earth’s  surface,  as  well 

i 

as  deriving  Snell's  law. 

Euler-Lagrange  Equations  of  a  Ray 

Consider  a  curve  C  joining  two  points  A  and  3  as  shown  in  Fig.  (l). 

The  time  t  required  for  a  wave  surface  of  constant  phase,  or  phase  front, 

to  travel  along  C  with  phase  velocity  v^  is 

B 

t  -  J'(l/v  )  ds  (18) 

A  P 


Fig,  1.  Variation  of  a  curve  C  between  two  points 
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where  ds  is  an  increment  of  axe  length.  Substituting  v^  =  c/n  from  2q. 
(l)  for  a  given  medium  into  3q.  (13)  gives 


3 

t  =  (l/c)  S  nds 

A 


(19) 


The  distance  that  the  phase  front  would  travel  in  the  medium  is  called 
the  optical  path  length  L,  and  is  given  as 


3 

L  =  ct  <=  S  nds 

A 


(20) 


Fermat’ s  principle  is  stated  as  follows!  >;hen  light  travels  from 
a  point  A  to  a  point  3,  it  travels  along  a  path  for  which  the  optical 
path  length  has  a  stationary  value.  From  the  calculus  of  variations 
this  "stationary  value"  is  usually  expressed  as  a  variation  6  of  the 
path  length  integral,  which  is  taken  to  be  a  minimum.  Thus  Fermat's 
principle  may  be  given  by 

3 

6  L  =  6  j  nds  «  0  (21) 

A 

Suppose  the  curve  C  joining  points  A  and  3  is  defined  by 

x(g),  y(g),  a(s)  (22) 

where  g  is  some  parameter  along  the  curve.  The  element  of  arc  length 
ds  along  C  is  then 


where  the  prime  denotes  differentiation  with  respect  to  g.  Substitut¬ 
ing  2q.  (23)  into  Eq.  (20 )  gives 


3 

L  =  J  wdg 
A 


(24) 


T 


where 


w  =  w(x,y,z,x' ,y' ,s* )  =  n(x,y,z)LX’2  +  y*2  +  z'2]  (2 l 

Consider  now  a  neighboring  curve  C  having  the  same  end  points  a 
and  3,  which  is  also  defined  in  terns  of  g  as  shown  in  Fig.  (l).  The 
variation  of  w  in  passing  from  a  point  P  on  the  curve  C  to  a  point  ?’ 
on  curve  C'  which  corresponds  to  the  same  value  of  g  is 


dw  ,  ,  \  dw  r  , 

=!  L 


where  the  notation  x^,  x2,  x^,  x£,  x^,  x^  is  used  to  denote  x,  y,  z,  x' , 
y' ,  z'  respectively.  Taking  the  variation  of  Eiq.  (24) ,  noting  that  the 
integral  sign  and  variational  symbol  6  are  commutable ,  and  substituting 
3q.  (26),  6L  becomes 

3  3 

6L  ■  6  j  ”dc  *  j  (&w)dg 
A  A 


A  Li=l  1  i-1  J 


Since  the  derivative  d  and  the  variation  6  are  likewise  commutable,  then 


6x^  «  dCdx^dg)  =  (d/dg)(6xi) 
Substituting  Eq.  (28)  into  Eq.  (27)  giv  es 


■J  L\*i*L*A 4^ 

A  Li=T  1  i=l  1  J 


Integrating  the  second  term  on  the  right  hand  side  of  the  equation  by 
parts  yields 
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;  3  3 

J  ^  lx^  df  ds  "  ^ 


A  1-1 


■I  X>«;>». 


A  1-1 

Since  the  curves  C  and  C*  have  common  end  points  A  and  3,  then  the  vari¬ 
ations  6k  ^  vanish  at  the  end  points.  Thus  the  first  term  on  the  right 
hand  side  of  Eq.  (30)  vanishes,  which  results  in  Eq.  (29)  being  written 
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A  1-1 

Fermat's  principle  requires  that  this  variation  must  be  zero  if 
curve  G  is  the  path  along  which  light  travels  from  point  A  to  point  B. 
That  is,  Eq.  (21)  must  be  satisfied  for  curve  C  to  be  a  ray.  Setting 
Eq.  (31)  equal  to  zero  and  reversing  the  terms  gives 


■\ t  [< 


A  i-1 


Since  all  were  introduced  as  completely  generalized  functions,  the 
Integrand  of  Eq.  (32)  must  vanish  to  avoid  trivial  solutions  in  which 


all  fix^  *  0  or  dg  -  0.  Thus  the  component  differential  equations 

dg  'Bx* '  dx  u 

d  fdv  \  dw  n 

dg  '3y* J  "  By 

(Hi  )  .  £"  „  0 

dg  'dz1 '  Zz  v 


(33a) 


(33b) 


(33c) 


IB 

must  be  satisfied  by  a  ray,  which  are  also  known  as  the  Suler-Lagrange 
equations . 

3qs.  (33)  are  valid  for  propagation  over  short  distances  where  the 
effects  of  earth  and  atmospheric  curvature  are  nil.  However  for  longer 
ranges  (>  100  nautical  miles)  it  is  necessary  to  consider  the  more  gen¬ 
eral  case  of  a  spherical  earth  and  a  spherically  stratified  atmosphere. 

Let  two  points  in  space  P  and  P'  be  located  a  short  distance  apart. 
The  arc  length  ds  between  the  points  is  given  in  spherical  coordinates  as 


Algorithm  Development 

Having  derived  the  general  equations  of  a  ray  in  a  spherically  strat¬ 
ified  medium,  it  is  now  desirable  to  solve  for  ray  height  versus  distance. 
For  the  sake  of  simplicity  in  computational  as  well  as  graphical  represent¬ 
ation  of  the  ray,  consider  any  arbitrary  0=  constant  plane  in  the  atmos¬ 


phere  as  shown  in  Fig.  (2).  Let  the  ray  lie  in  the  constant  plane  and 
let  the  refractive  index  vary  as  a  function  of  radial  distance  r  for  r  > 
a,  where  a  is  the  earth's  radius. 


Fig.  2.  Bay  geonetry  in  polar  ooordinatoa. 


Substituting  g  -  r  into  Bqs.  (3*0  and  (35)  givoe 
[l*<rg)2]1/2«r 

.  -  «<r)  [l  ♦  (r  g)2]  1/2 
Inserting  Bq.  (38)  into  £qs.  (36)  gives 


ftros  the  geoastry  of  Fig.  (2)  it  should  be  noted  that 

r  -  r(e) 
e  -  e(r) 

However,  Sq.  (38)  is  an  explicit  Amotion  only  of  r  and  d8/dr, 
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vfcere  the  prise  denotes  differentiation  with  respect  to  r  (or  g).  The 
derivative  of  w  is  written 

dw  -  f|  dr  +  |*  d0'  (42) 

Thus  while  the  total  derivative  dw/de  say  exist  for  Eq.  (42),  the  partial 
derivative  dw/dd  does  not.  Therefore 

ft  -  0  «» 

Combining  Eqs.  (39)  and  (43)  yields  an  Buler-Lagrange  equation  in  polar 


coordinates  given  by 


d 

dw 

3r 

k>J 

m  0 


Similarly  letting  g  -  0  in  Eqs.  (34)  and  (35)  gives 


ds 


*[<Si>2+'2] 


1/2 


- 

I 


-  -  »<r)  I  (f|)  +  *2| 


de 

,11/2 


J 


Substituting  Eq.  (46)  into  Eqs.  (36)  yields 


-g-° 


d 

dw 

de 

(44) 


(45) 

(46) 

(47) 


In  this  case  Eq.  (46)  is  an  explicit  function  only  of  r  and  dr/de, 


or 

w  -  w(r,  r’ )  (48) 

where  the  prise  indicates  differentiation  with  respect  to  0  (or  g).  The 
derivative  then  becomes 

dw  -  fj  dr  +  dr'  (49) 

Since  both  dw/dr  and  dw/dr  exist  far  Eq.  (49),  then  Eq.  (47)  remains  un¬ 
changed,  thus  providing  a  second  Buler-Lagrange  equation. 
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Eq.  (44)  is  useful  since  it  produces  a  form  of  Snell' s  law  for 
spherically  concentric  refractive  media,  and  Eq.  (47 )  is  important  be¬ 
cause  it  results  in  a  second  order  differential  equation  which  defines 
the  trajectory  of  a  ray  in  a  polar  coordinate  system. 

Substituting  Eq.  (38)  into  Eq.  (44)  and  integrating  Eq.  (44)  gives 


»<g> 


n(r)Ll  +  (r 


2  1/2" 

— )  1  =  K. 

dr '  -I 


(50; 


which  becomes 


Aw  s 


& +  c*  in 


2  1/2 


=  K 


(51) 


where  K  is  a  constant  of  integration.  Prom  Fig.  (2)  it  may  be  seen  that 


sin  £  = 


rde 


[(dr)2  +  (rde)2]1^2 
_ 


dr 


&  ♦  (r  f )  ] 


TTTz 


(52) 


where  £  is  the  angle  between  the  ray  and  the  radius  vector  r.  Substitut¬ 
ing  Eq.  (52)  into  Eq.  (51)  gives 

m(r)  sin  £  «  K  (53) 

which  is  Snell's  law  for  spherically  concentric  media.  Evaluating  K  at 
some  known  point  (e.g.,  the  emitter)  yields 

rn(r)  sin  =  ron(rQ)  sin  (54) 

where  the  zero  subscript  refers  to  known  values  at  the  emitter. 

Substituting  Eq.  (54)  into  Eq.  (52),  solving  for  de/dr,  and  inte¬ 


grating  gives 


Sq.  (55)  could  be  used  to  calculate  the  ray  trajectory  except  for  the 
fact  that,  as  Fig.  (3)  illustrates,  0  may  become  a  double  valued  func¬ 
tion  of  r. 


TO  CENTER 
OF  THE  ERR.TH 

Fig.  3.  Ray  trajectory  for  a  double  valued  solution  to  Eq.  (55). 

An  initially  rising  ray  is  refracted  downward  toward  the  earth  resulting 
in  two  possible  values  of  0  (0^  and  0^)  for  a  single  value  of  r  =  r^. 

A  more  satisfactory  approach  to  determining  the  ray  trajectory  is  to 
expand  Eq,  (4?)  since  it  permits  the  solution  of  r  as  a  function  of  0. 
Substituting  3q.  (46)  into  Eq.  (47)  gives 


-  -  r2]1/2l-  0 


(56) 


Carrying  out  the  differentiation  and  rearranging  the  terms  yields 


a2r  _  2/drs2  ,  r2  aCn(r)]  f/1  dr 
2  ”  r^dG'  n(r)  Sr  [_^r  dG 


)  +  1  +  r 


Substituting  the  relationship 


a?  &»(«>:  -  \  Tr 


into  So.  (57)  gives 


^ - ? (i)2 + *2 &*&<*>:»  r^i)2-] +  -  <»> 

dG  L  J 

Eq.  (59)  is  more  useful  when  it  is  changed  from  polar  coordinates 
(r,  G)  to  the  coordinates  of  distance  over  the  earth's  surface  and  alti¬ 
tude  above  the  earth's  surface  (x,  h).  Letting 


r  =  a  +  h 


e  =  * 


(60a) 

(60b) 


Eq.  (59)  may  be  written 


d  h  _  (a  +  h 


^  [  c^h)  g]: 

mm 


+  (a  +  h)  Ih  (lnCn(h)3) 


([V^s]  +1]  +1< 


Eq.  (6l)  forms  the  basis  for  all  ray  trajectory  computations  through¬ 
out  this  paper,  and  is  in  agreement  with  the  ray  equations  developed  by 
22 

.-artree  et  al  for  use  on  an  analog  differential  analyzer.  The  results 
presented  in  this  thesis  are  obtained  by  numerically  integrating  Eq.  (6l) 
on  a  digital  computer  by  using  a  standard  fourth  order  Runge-Kutta  inte¬ 
gration  algorithm. 


I5ow  that  the  ray  trajectory  may  be  found,  it  is  useful  to  derive  a 
solution  for  the  time  of  propagation  along  the  ray  path.  Eq.  (19)  gives 
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the  tiae  it  takes  for  a  wave  flront  to  travel  between  two  points  along  a 
specified  curve  C.  If  C  is  the  ray  path  which  satisfies  the  Euler- 
Lagrange  equations  in  a  polar  coordinate  system  (Eqs.  (36)  for  a  <£  - 
constant  plane)  then  substituting  Eqs.  (37)  and  (45)  into  Eq.  (19)  yields 


*c/“(r)[ 


1  +  <*f!) 


211/2 


S  [(g)2  *  r2]  lA  d« 


where  (rQ,  6q)  are  the  esitter  coordinates  and  (r,  e)  are  the  coordin¬ 
ates  of  a  point  on  the  ray.  Both  solutions  axe  valid  although  as  in  the 
case  of  Eq,  (55) »  Bq.  (62)  nay  be  double  valued  for  sons  ray  trajectories. 

Differentiating  Eq.  (63)  and  substituting  the  coordinate  relation¬ 
ships  of  Eqs.  (60)  gives 

dt  „  1  ,  .  f  ,dhv2  ,  /a+hx2!1/2 

dx  c  a'h'  ^dx^  +  ^  a  )  J  (64) 


Bq.  (64)  is  useful  since  it  provides  the  phase  iufaraation  necessary  to 
calculate  the  total  field  in  regions  of  signal  interference.  It  is  In¬ 
tegrated  with  Eq.  (6l)  to  provide  the  position  and  tine  coordinates  (x, 
h,  t)  of  each  ray. 


CHAPTER  IV 

PROPERTIES  OF  GEOMETRIC  OPTICS 


'.'hile  the  method  of  geometric  optics  is  a  simple  technique  for 

solving  many  propagation  problems,  there  exist  a  number  of  conditions 

which  limit  its  useage.  One  way  of  presenting  these  restrictions  and 

other  related  properties  is  to  examine  the  wave  equation  and  its  solu- 

1  1 5 

tion  in  both  homogeneous  and  inhomogeneous  media.  * 

Cave  and  Tikonal  Equations 

Consider  a  nearly  perfect  dielectric  medium  that  is  charge-free  (no 
sources  exist),  unbounded  (infinite  in  extent),  linear  (the  permeability 
p.  and  permittivity  €  of  the  medium  are  independent  of  the  fields),  homo¬ 
geneous  (u.  and  €  are  independent  of  position),  and  isotropic  (p.  and  € 
are  scalar).  The  wave  equation  is  then  given  as 

V2  3  +  k2  E  -  0  (65) 

where  E  is  the  electric  field  intensity  vector  and  k  is  the  propagation 
constant  of  the  medium.  Prom  elementary  electromagnetic  theory  it  is 
known  that  for  a  nearly  perfect  dielectric  medium 

k  «  2v/\  “  CO/v_  ■  CO  "\f  d-  (66) 

where  \  is  the  wavelength,  CO  is  the  radian  frequency,  and  v^  is  the 
phase  velocity  of  the  field  in  the  medium.  The  free  space  propagation 
constant  is 

ko  *  2Tr/\>  "  00 /c  "  CO  V^o-o  (6?) 
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vhere  the  zero  subscript  denotes  free  space  values  and  c  is  the  speed  of 
light  in  vacuum. 

Substituting  3qs.  (2)  and  (3)  for  a  non-magnetic  medium  into  Eq. 


(66)  gives 


'WV^oVr  =  kc 


where  €  is  the  relative  permittivity.  The  wave  equation  may  now  be 


written  as 


r2  E  +  k? 


k~  €  2  «  0 
o  r 


Finally,  substituting  So.  (4)  into  Eq.  (63)  yields 


k  *  k  n 
o 


and  Eq.  (65)  becomes 

V2  E  +  k2n2  2-0  (71) 

Thus  3qs.  (65),  (69),  and  (71)  give  three  forms  of  the  wave  equation 
for  2  in  terms  of  the  propagation  constant,  relative  permittivity,  and 
index  of  refraction  of  the  medium. 

So.  (71)  is  particularly  interesting  since  it  explicitly  contains 
the  refractive  index.  Using  rectangular  coordinates  to  simplify  the 
treatment  of  plane  waves,  Eq.  (71)  may  be  expressed  in  terms  of  its 
scalar  comoonent  fields 


V2  Ex  +  k2  n2  2x  -  0 


(72a) 


T2  Ey  +  k2  n2  Ey  -  0 
V2  Sz  +  k2  n2  Ez  -  0 


(72b) 


(72c) 


?or  a  homogeneous  medium  where  n  is  constant,  the  plane  wave  solution 
for  Sqs.  (72)  will  have  the  form 
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i 


j 

) 


i 


E  -  *.'J(k!»n)'r 


(73) 


Ae 

where  kQ  Is  the  free  apace  propagation  vector  which  is  noraal  to  the 

wavefront  and  pointed  In  the  direction  of  notion  of  the  planes  of  con* 

A  ■*  * 

stant  phase,  n  is  the  unit  vector  in  the  direction  of  kQ,  and  r  is  the 

position  vector  of  the  point  in  space  at  which  the  field  is  to  be  cal¬ 
culated. 


In  reality,  atmospheric  refractive  index  is  not  a  constant  but 
varies  slightly  with  position  so  that  Eq.  (73)  is  no  longer  valid.  How¬ 
ever  a  solution  to  the  wave  equation  does  exist  in  a  farm  very  similar 
to  Eq.  (73).  Replacing  A  and  n  by  Q  and  S,  which  are  real  functions  of 
position,  the  solution  be cones 

E  -  Qe^o3  (74) 

Substituting  Eq.  (?4)  into  any  of  the  component  fields  of  Eqs.  (72) 
yields 

[v2^  -  k2Q(VS)2  -jk^S  -  j2k0  VQ»VS>_jko3  +  k2nZQe'^oS  -  0  (75) 

Rearranging  the  terns,  dividing  by  Qk02,  and  dropping  the  exponent,  this 
be cones 


n2  -  (VS)2  -  (j/k^E^S  ♦  2VQ-VS/q3  -  vfyQk02  (76) 

Mow  if  kQ  is  sufficiently  large  such  that  the  following  two  conditions 
are  satisfied i 


v« 

h  o 


(77) 


'j 

i 


(78) 
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then  Bq.  (?6)  reduces  to 

n2  -  (VS)2  (79) 

which  is  called  the  equation  of  the  eikonal.  While  Eq.  (79)  nay  be  used 
directly  to  construct  the  ray  paths  of  geometric  optics , ^  it  will  be 
used  here  to  derive  the  relationship  of  ray  density  to  the  power  density 
of  the  field  and  to  present  the  limitations  of  ray  optics. 


The  Jacobian  and  Power  Density 

Consider  the  transformation  of  a  point  in  space  between  the  coor¬ 
dinates  (x,  y,  z)  and  (u(  v,  w).  Let  the  coordinates  (x,  y,  z)  be  de¬ 
finable  in  terms  of  (u,  v,  w)  such  that 

x  -  x(u,  v,  w)  (80a) 

y  -  y(u,  v,  w)  (80b) 

s  -  a  (u,  v,  w)  (80c) 

If  at  a  given  point  PQ  where  (xq,  yQ,  zQ)  -  (uq,  vq,  wq)  the  Jacobian 


of  the  transformation 


d(u,  v,  wj 


dx  dx  dx 
du  dv  dw 

*1  *1  *1 
du  dv  dw 

dz  dz  ds 
du  dv  dw 


is  not  zero,  and  if  all  the  partial  derivatives  in  Bq.  (81)  are  contin¬ 
uous  throughout  some  region  which  includes  PQ,  then  Sqs.  (80)  may  be 
solved  uniquely  for  u,  v,  and  w  in  terms  of  x,  y,  and  z  in  this  region. 
Holding  each  of  the  variables  u,  v,  and  w  constant  in  turn  gives  three 
parametric,  or  coordinate,  surfaces  through  a  point  in  space,  further¬ 
more,  holding  any  pair  of  the  variables  constant  yields  a  parametric,  or 
coordinate,  curve  pasting  through  the  point. ^ 


Consider  next  a  volume  element  dY  bounded  by  a  bundle  or  "pencil"  of 


rays  and  surfaces  of  constant  phase  as  shown  in  Fig.  (4).  Let  the  vol¬ 
ume  he  defined  by  the  surfaces  u  -  ,  u  *  +  du,  v  =  v.,  ,  v  =  v^  +  dv, 

v  =  w^(  and  w  =  w^  +  dw,  where  the  variables  u,  v,  and  w  are  solvable 
within  the  volume  such  that 

u  -  u(x,  y,  z)  (82a) 


v  *  v(x,  y,  z) 
w  =  w(x,  y,  2) 

If  the  surfaces  of  constant  phase  are  given  by 


(82a) 

(82b) 

(82c) 


w  »  w„ 


w  =  w^  +  dw  «=  +  dS  (84) 

where  3^  and  dS  represent  the  phase  angle  and  phase  angle  increment  for 
3q.  (74),  then  Eqs.  (82b)  and  (82c)  produce  two  families  of  surfaces 
whose  intersections  are  the  rays  between  S^and  dS.  It  should  be  noted 
that  these  families  are  not  uniquely  determined  by  the  ray  pattern  since 
it  is  possible  to  construct  an  infinite  number  of  families  having  the 
same  intersections  as  Eqs,  (82b)  and  (82c). 

Equating  the  imaginary  part  of  Eq.  (76)  to  zero  yields 


Combining  Eq.  (85)  with  the  vector  relationship 


v.(q2vs)  «  2Q(VQ.ys)  +  q^s 


gives 


V.(tVS)  -  0  (87, 

2 

Since  Q  73  is  solenoidal  (i.e.,  its  divergence  is  zero),  there  exists  a 
vector  U  such  that 


Q2VS  -  V  x  u 


(88) 
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Fig.  4.  The  volume  element  dy  bounded  by  rays  and  surfaces  of  con¬ 
stant  phase.  The  face  abed  lies  on  the  surface  w  ■  S-,and  the  face  efgh 
lies  on  the  surface  w  »  S.  +  dS.  Similarly  abfe  and  aegh  are  portions  of 
the  surfaces  v  -  v.,  and  v  -  v.  +  dv,  respectively.  The  remaining  two 
faces  adhe  and  begf  lie  on  the1  surfaces  u  *  u^,  and  u  -  +  du,  respect¬ 

ively.1 

where  the  curl  of  U  lies  in  the  direction  of  VS.  By  correctly  choosing 
U,  the  surfaces  u  and  v  may  be  found  which  satisfy  the  wave  equation  and 
the  rays  within  d  y,  although  surface  families  other  than  those  of  u  and 
v  nay  also  exist. 

It  was  stated  in  Chapter  I  that  rays  are  perpendicular  to  the  wave- 

fronts,  or  surfaces  of  constant  phase,  which  are  propagating  in  an  iso- 
16 

tropic  medium.  Since  V3  defines  the  normal  to  a  wavefront,  then  for 
the  surfaces  of  u  and  v  to  intersect  along  the  ray  paths,  Vu  and  Vv  must 
be  normal  to  V3,  and  the  vector  Tux  W  must  be  parallel  to  VS.  This  *' 
latter  vector  is  used  in  the  identity 

V  x  (uVv)  ■  Vu  x  W  +  u£v  x  (Vv)3 

-  Vu  x  W  (89) 

where  the  term  u£v  x  (Vv)J  vanishes  since,  from  vector  calculus,  the  curl 
of  a  gradient  is  zero  when  continuous  mixed  partial  derivatives  are  as¬ 
sumed  so  that  the  order  of  differentiation  is  immaterial.  Allowing  u  and 
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v  to  be  veil  behaved  functions  with  continuous  derivatives,  Sq.  (89) 
rives 

U  -  uVv  (90) 
where  tne  curl  of  U  lies  in  the  direction  of  VS.  Substituting  Eq.  (90) 
into  Sq.  (88)  yields 


3^3  *  V  x  (uW)  ■  Vu  x  W  (91). 

Tailing  the  dot  product  of  each  side  of  Sq.  (91)  with  VS  and  substitut¬ 
ing  Eq.  (79)  gives 


,2  m  (Vu  x  Vv).V3 


for  which  the  numerator  may  be  written 


du  5u  &u 
dx  dy  oz 


(92) 


(Vu  x  Vv). V3 


dv  dv  dv 
dx  dy  dz 

dw  dw  dw 
dx  dy  dz 


(■u.v«w)  •  I 

(x.y.z)  J 


/no') 

\7J) 


where  J  is  the  Jacobian  of  x,  y,  and  z  with  respect  to  u,  v,  and  w. 

Thus  Sq.  (92)  may  be  written 

Q2  -  ±2  (94) 

The  volume  increment  of  dy  may  be  expressed  in  terms  of  the  Jacob-  „ 
lan  by 


dy  -  J(u,  v,  w)  dudvdS  (95) 
let  dA  be  the  element  of  area  on  the  surface  w  «  51  which  is  bounded  by 
the  surfaces  u  «  constant  and  v  -  constant  as  shown  in  Fig.  (4).  Since 
tne  rays  are  normal  to  their  respective  wavefronts,  Sq.  (95)  may  be  ap¬ 
proximated  by 


B  ■  (i  H  (lOrb) 
For  a  plane  sinusoidal  wave  propagating  through  an  infinite  homogeneous 
or  inhomogeneous  medium,  3a  c.  (100 )  dictate  that  S  and  H  are  transverse 
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or  perpendicular  to  the  propagation  vector  k  of  the  field,  and  normal 
to  each  other.  fVom  Eos.  (100)  and  (101 )  this  relationship  is  written 


€/u  ft  x  3 


(102) 


where  ft  is  the  unit  vector  in  the  direction  of  k. 


The  density  of  the  outward  flux  of  energy,  or  power  density,  far  a 
time-varying  field  passing  through  a  closed  surface  is  given  by  the 
Poynting  vector  ?,  where 


P  «  I  x  H  (103) 

Following  the  standard  convention  that  physically  meaningful  electric 
and  magnetic  fields  are  represented  by  the  real  part  of  complex  quanti- 
ties  (i.e. ,  cos©  *  He[eJ  j),  the  time-averaged  Poynting  vector  for  sin¬ 
usoidal  fields  is  written 

“  2  "  2  X  **-  '  (104) 
where  H*  denotes  the  complex  conjugate  of  H  and  the  operator  3e  means 

"the  real  part"  of  the  complex  Poynting  vector  F  *  3  x  H*.  Substitut- 

v 

ing  Eq.  (102)  into  Eq.  (104)  gives 

rav  -  I  V5V  I S I  2  ft  (105) 

which  may  be  integrated  to  determine  the  average  outward  energy  flux  or 
average  power  p  through  a  closed  surface  S’ ,  such  that 


(106) 


o 

Substituting  Sq.  (99)  into  Eq.  (94)  and  using  n  -  VS* VS  from  Eq. 
(79)  yields 


dudv  «  Q2ndA  -  Q2V3-(VS/n)  d A  »  £2VS*ftdA  (107) 

where  V3/n  ■  n  is  the  unit  vector  normal  to  the  wavefront  at  the  surface 
element  dA.  Integrating  Eq.  (107)  over  a  closed  surface  S'  which  includes 
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the  wavefront  surface  3^  gives 

ijdudv  *  $$  Q^VS  *  ndA 
S'  S'. 

=  ®  Q27S.d^  (108) 

S’ 

2  •*  2 

Noting  from  Sq.  (?4)  that)Q  VSjis  proportional  to  I  El  and  lies  in  the 
direction  of  n,  then  a  comparison  of  Eqs.  (105),  (106),  and  (108)  indi¬ 
cates  that^Jdudv  is  the  total  outward  power  radiated  through  the  surface 

S 1 

S' .  Thus  the  element  dudv,  which  is  "bounded  by  the  ray  bundle  of  Fig. 
(4),  is  the  time-averaged  energy  flux  of  the  Poynting  vector  P  *  Q  V3 

ctv 

through  the  wavefront  surface  increment  cLA,  or 

dudv  -  Q2V5 •  ndA  -  P  *ndA  (109) 

It  is  this  relationship  of  the  ray  bundle  element  dudv  to  the  average 
power  density  at  the  wavefront  which  is  used  in  the  power  density  and 
field  strength  calculations  presented  in  Chapter  VI  of  the  thesis. 


Restrictions  of  Ray  Theory 

The  restrictions  in  the  use  of  geometric  optics  result  from  the  in¬ 
equality  of  Eq.  (77).  Using  the  vector  identity 


p.(23)  .  +  VQ-V(~) 

vnQ'  nQ  H  nnQ' 


(no) 


Eq.  (77)  may  be  written 

-A  „  _JL  v.(_SLu_l_  (S).(-2SL)  +  (  _22_).(_22_)  «  1  /vn) 

1 2„20  k„n  v  nQ'  k  n  '  n'  'k  nQ'  '  k  nQ'  'k  nQ'  1 
00  o  ~ 


kJTQ  -o- 


-0n®  V5' 


Taking  the  square  root  of  Eq.  (94)  and  substituting  it  into  VQ/k  nQ  gives 


-22 _ 


k0n2 


V(J1/>2)  ,  1  Vn 


Lvc^2) 


k  n  n 


(112) 


Combining  Eqs.  (Ill)  and  (112)  yields 
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k*n2Q 


I 

* 

> 

JL(ffi) 

k  rr  nJ 

0 

4, 

1 

<1 

> 

0 

fo 

v — . 

1 _ 

e 

1 

N»1 

1 _ 

kon 


Lko 


ViiH L 

nJW2) 


kon 


JL/2&) 

k  n'  n' 
o 


«  1 


(113; 


which  is  satisfied  only  when 
X  IVhl 


V 


«i 


(114) 


and 


«  1 


(115) 


1  kj1/2)! 

V  (W2) 

Furthermore,  Eq.  (H3)  requires  that  (l/kQn)(Pn/n)  and  (l/kon)Cv(J^^2)/ 
be  smooth  so  that  their  divergences  are  likewise  small* 

Eqs.  (114)  and  (115)  specify  the  conditions  governing  the  use  of 
ray  optics.  Recalling  that  kQ  »  2Tr/XQ  and  n  *  kQ/k  -  X0/l,  Eq.  (114)  may 
be  written 


°-(^)  -  -£(^i)  «  1 


(116) 


2rmv  n'  2rr  n' 

which  states  that  the  relative  change  in  the  refractive  index  over  a  dis¬ 
tance  approximately  equal  to  a  wavelength  of  radiation  must  be  small  com¬ 
pered  with  unity.  Thus  Eq.  (114)  becomes  a  better  approximation  as  wave¬ 
length  decreases,  and  is  perfectly  satisfied  in  a  homogeneous  medium 
where  7a  •  0,  On  the  other  hand,  a  medium  with  sharp  discontinuities  in 
the  refractive  index  violates  Eq.  (114)  at  the  discontinuities.  Therefore 
while  rays  obey  Snell's  law  on  either  side  of  a  discontinuity  in  the  re¬ 
fractive  index  (e.g. ,  either  side  of  the  boundary  between  two  dielectric 
ftedla),  the  region  at  the  discontinuity  itself  cannot  be  described  by 
geometric  optics,  but  Instead  must  be  treated  by  the  more  general  wave 


i 

i  or  physical  optics  theory. 

i 

i 

3q.  (114)  also  implies  that  the  radius  of  curvature  of  a  ray  oust 
be  much  greater  than  a  wavelength  of  radiation.  Consider  a  ray  passing 
from  point  A  to  point  A'  as  shown  in  Fig.  (5).  Let  the  refractive  in¬ 
dex  vary  as  a  function  of  radial  distance  r  from  the  center  of  the 
earth,  and  let  p  be  the  radius  of  ray  curvature  between  A  and  A' . 

Uext  consider  a  snail  wavefront  increment  A3  which  is  normal  to  the  ray 

at  noint  A  at  time  t  .  The  wavefront  arrives  at  A’ 3'  when  the  time  is 

0 

t  +  t.  Since  the  wavefront  must  travel  from  A  to  A'  in  the  same  anoun-l 
o 

of  time  t  that  it  takes  to  travel  from  3  to  3* ,  the  arc  lengths  between 
the  two  paths  nay  be  written  as 


vt 


(117; 


and 

C  +  d<J  *=  (v  +  dv)t 

where  v  is  the  phase  velocity  between  A  and  A* . 
tween  the  two  wavefront  positions  becomes 

\lt  .  0  ,  *g 
T  p  p  +  dp 

v£_  (v  +  avH 

■p-p.ap 


(us; 

Thus  the  angle  \jj  he¬ 


ar  slml, 


Substituting 


from  3q.  (l)  and 


djD  dv 
P  "  v 


c 

V  C  — 

n 


(119) 


(1207 


(121) 


dv  *  -  —  dv 
n^ 


(122) 
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Tig.  5.  leone try  for  ray  curvature  in  polar  coordinates 
into  iq.  (120)  jives 

f--~>  <*» 

Iron  ”i§*  (5)  it  may  be  seen  that  the  incremental  radial  displacement 
dr  of  the  ray  as  measured  from  the  center  of  the  earth  is  approximately 

dr  •  dp  •  coso:  (12*0 

"hers  0;  is  the  angle  between  the  earth  radius  vector  r  and  the  wave- 
front  .  c.  is  also  the  angle  of  ray  elevation  with  respect  to  the  local 
horizontal  plane  (i.e,,  the  plane  normal  to  r)  drawn  at  point  A,  Solv¬ 
ing  for  dp  in  dq.  (124)  and  substituting  into  So.  (123)  yields 


cosa 


(125) 


~ir.ce  the  refractive  index  is  a  function  only  of  r,  dr. /dr  may  be  re¬ 


placed  by  rn  so  that 


^  =  — —  cos  a 

P 


Thus  the  requirement  ox  Tq.  (114)  that  Tn/n  be  small  as  compared  to  a 
wavelength  of  distance  forces  the  radius  of  ray  curvature  in  To.  (126) 
to  be  large  over  the  same  distance.  Furthermore,  Tq.  (126)  clearly  in¬ 
dicates  that  ray  bending  is  dependent  upon  the  relative  change  in  the 
refractive  index  of  the  medium.  For  a  ray  to  follow  tne  curvature  of 
the  earth  as  in  the  case  of  wave  trapping  or  ducting,  To.  (12c)  yields 


(127) 


where  p  =  a  is  the  earth's  radius,  n  *  1  for  air,  and  a  =  0.  In  free 


space  where  Vn  =  0,  then 


which  enables  the  rays  to  travel  in  straight  lines. 


f  hf  ' 


returning  to  the  second  major  condition  of  geometric  optics,  To. 

(99)  indicates  that  Z  is  directly  proportional  to  the  cross  section  cu> 

i  /? 

of  the  ray  bundle,  thus  implying  that  is  related  to  the  spacing  be¬ 

tween  neighboring  rays.  Therefore  To.  (115)  states  that  the  fractional 
change  in  the  spacing  between  neighboring  rays  over  a  wavelength  of  dis¬ 
tance  must  also  be  small  compared  with  unity.  This  condition  is  violated 
where  rays  either  converge  or  diverge  rapidly  as  is  often  the  case  along 
the  edge  of  a  reflected  or  refracted  envelope  of  rays.  Tuch  a  region, 
known  as  a  caustic,  is  generally  characterised  by  a  ray  density  which  in¬ 
creases  steadily  up  to  the  edge  of  the  ray  envelope  and  then  sharply  de¬ 
creases  to  zero  beyond  the  envelope  boundary.  Furthermore,  Tc.  (115)  is 
violated  whenever  the  cross  section  of  a  ray  pencil  vanishes  as  it  does 
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when  the  pencil  emerges  from  a  point  source  or  passes  through  a  focal 
point  or  caustic.  Thus  point  source  regions,  focal  points,  arc!  caustics 
must  be  excluded  from  any  meaningful  consideration  in  ray  optics,  al¬ 
though  like  discontinuities  in  the  medium,  they  may  be  handled  by  wave 
theory.1*  1^»  1*^ 

Finally,  since  ray  curvature  is  affected  only  by  the  relative  change 
in  the  refractive  index  and  the  angle  of  ray  elevation  in  the  medium, 
there  is  no  mechanism  within  classical  geometric  optics  to  account  far 
the  phenomenon  of  wavefront  diffraction  around  the  edges,  earners,  or 
vertices  of  boundary  surfaces.  Modifications  to  the  theory,  which  are 
generally  referred  to  as  the  Geometrical  Theory  of  Diffraction  (GTD), 
have  been  made  in  recent  years  to  overcome  this  deficiency.  ^  While  the 
GTD  lies  beyond  the  scope  of  the  thesis,  it  is  worth  noting  that  this 
extension  of  ordinary  ray  optics  is  also  based  upon  Fermat's  principle 
and  the  equation  of  the  elkonal  in  modified  farm. 


CHAPTER  V 


COMPUTER  SIMULATION  DESCRIPTION 

General  Characteristics 

The  previous  discussion  of  refractivity  and  geometric  optics  the¬ 
ory  provides  a  model  of  radio  wave  propagation  through  a  non- ionized 
and  time  invariant  atmospheric  medium  that  is  radially  stratified  with 
respect  to  the  earth's  center.  The  model  follows  the  polar  coordinate 
convention  of  Fig.  (2),  which  assumes  azimuthal  symmetry  in  the  atmos¬ 
phere  about  the  emitter.  A  computer  program  which  uses  this  model  was 
developed  for  the  thesis  to  serve  as  an  analytic  tool  in  comparing  the 
field  calculated  from  ray  theory,  with  available  experimental  data  for 
an  isotropic  emitter.  A  description  of  the  program  and  its  algorithms 
is  presented  here,  with  specific  details  regarding  program  useage  given 
in  the  Appendices, 

Program  ATREF  (atmospheric  refractivity)  is  written  in  Extended 
FORTRAN  language  for  the  CEC  5600  digital  computer  system.  It  requires 
approximately  65,000  octal  words  of  memory  and  15  seconds  of  central  pro 
cessor  time  for  execution  of  a  run.  Output  is  provided  by  a  standard 
132  column  line  printer  and  an  on-line  GALGOMP  plotter. 

Given  the  following  input  in  card  image  formi 

a)  Atmospheric  refractivity  versus  altitude 

b)  Emitter  parameters  (frequency,  polarization,  pulse  width,  beam 
upper  and  lower  angle  limits,  altitude) 

c)  Earth  surface  parameters  (land  or  sea,  surface  height  variation 
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the  program  calculates  the  ray  trajectories  and  the  time  of  propagation 
along  each  ray  path  from  Eqs.  (6l)  and  (64),  respectively.  The  rays  are 
launched  in  the  vertical  (i.e. ,  elevation)  plane  with  a  0.02  degree  angu¬ 
lar  separation  to  ensure  a  reasonable  ray  density  at  long  ranges.  The 
program  allows  the  rays  to  be  launched  within  a  beam  whose  total  angular 
width  is  a  maximum  of  one  degree,  although  this  beam  width  may  be  enlarged 
by  increasing  the  array  dimensions  and  thus  the  computer  memory  require¬ 
ment. 

The  coordinates  of  distance  along  the  earth's  surface,  height  above 
the  surface,  and  time  of  propagation  (x,  h,  t)  are  stored  for  each  ray  at 
regular  intervals  of  dx.  These  coordinates  are  then  employed  in  an  algo¬ 
rithm  based  upon  Eq.  (109)  which  computes  the  given  emitter's  relative 
power  density  and  field  strength,  normalized  to  free  space  (vacuum)  values. 
Each  ray  that  intersects  the  earth  is  reflected  in  the  specular  direction 
and  is  assigned  a  complex  scattering  coefficient  which  is  a  function  of 
the  incidence  angle,  emitter  frequency  and  polarization,  terrain  type,  and 
surface  roughness.  The  scattering  coefficients  of  each  ray  are  then  en¬ 
tered  into  the  power  density  and  field  calculations  to  model  the  attenua¬ 
tion  and  phase  change  of  earth  reflected  signals. 

The  computer  simulation  provides  the  user  the  option  to  select  any  or 
all  of  the  following  items  for  output  in  either  printed  or  plotted  form, 
or  both i 

a)  Profile  of  refractivity  N  versus  altitude 

b)  Profile  of  the  vertical  gradient  of  N  versus  altitude 

c)  Ray  trajectories  (plotted  only) 

d)  Relative  field  strength  versus  altitude  and  distance  along  the 
earth's  surface 

e)  Relative  power  density  versus  altitude  and  distance  along  the 
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earth's  surface. 

The  plotted  output  of  options  (d)  and  (e)  shows  the  relative  field 
strength  (or  power  density)  of  the  total  E  field  in  space,  whereas  the 
printed  output  permits  a  more  detailed  breakout  of  the  total  field  and 
each  of  its  components,  as  in  the  case  of  multipath  interference  where 
the  field  in  space  is  the  resultant  of  the  fields  from  the  main  propa¬ 
gating  wave  plus  a  wave  component  which  has  been  reflected  from  the 
earth.  The  program  also  prints  the  time  difference  of  arrival  between 
the  main  and  multipath  wavefronts  (which  is  useful  in  some  radar  prob¬ 
lems)  and  the  number  of  wave  components  in  the  field  strength  arid  power 
density  calculations. 


Refractivlty  Models 

Uhile  the  simulation  is  generally  used  with  card  input  refractivity 
profiles,  the  user  may  select  either  of  two  reference  refractivity  mod¬ 
els  which  are  stored  in  the  simulation.  Both  are  based  on  the  Central 
Radio  Propagation  Laboratory  (CKPL)  model  given  by  Eq.  (l6).  Since  all 
calculations  in  the  simulation  are  in  the  units  of  feet,  the  constants 
in  Eq.  (l 6)  become 


Ng  -  313.0 

(129a) 

c  -  0.0000438  ft"1 
© 

(129b) 

for  the  stored  exponential  refractivity  model  and 

N  -  0.0 
s 

(130a) 

ce  -  0.0  ft"1 

(130b) 

for  the  stored  free  space  model.  If  the  program  is  used  with  an  input 
refractivity  profile,  a  curve  fit  is  made  to  Eqs.  (l6)  and  (129)  to  model 
atmospheric  refractivity  for  altitudes  which  lie  above  the  highest  point 


is  the  profile. 
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Bay  Algorithms 

The  simulation  uses  a  standard  fourth  order  Runge-Kutta  Integration 
algorithm  to  solve  the  ray  trajectory  and  time  of  propagation  equations 
given  in  Chapter  HI  by  Eqs.  (6l)  and  (64),  respectively.  The  integrat¬ 
ion  step  size  used  in  the  Runge-Kutta  algorithm  is  equal  to  one  tenth 
the  distance  along  the  earth's  surface  that  the  user  has  specified  far 
printing  or  plotting  output  data.  For  example,  if  the  user  wants  data 
printed  or  plotted  every  X  nautical  miles,  then  the  step  size  becomes 
X/lO  nautical  miles. 

Eqs.  (6l)  and  (64)  arc  written  in  state  variable  notation  for  num¬ 
erical  integration.  Letting  the  ray  altitude  above  mean  sea  level  be 
given  in  state  notation,  then  h  -  h^,  dh/dx  -  h^,  and  d2h/dx2  -  h^"  » 
h£.  The  second  order  differential  equation  of  Eq.  (6l)  becomes  a  pair 
of  first  order  state  equations  which  are  given  by 


a  + 


- 


2ftr+h^ 


(» +  [t(r^)  i-iD2  +  ij 2  + 1 


(131) 


and 


^  \  (132) 

where  n(hj)  -  N(hj)»l<T^  ♦  1  is  described  by  either  a  piece-wise  linear 
function  connecting  the  input  refractivity  profile  data  points,  or  the 
stored  function  of  Eq.  (l6).  Similarly,  Eq.  (64)  becomes 


(133) 


where  the  values  of  and  h^  are  obtained  from  the  integration  of  Sqs. 
(131)  and  (132). 


Power  Density  and  Field  Strength  Algorithms 
The  calculation  of  the  power  density  of  a  field  is  based  upon  Eq. 
(109)  which  states  that  the  power  density  is  directly  related  to  the 
element  of  flux  dudv  bounded  by  a  ray  bundle  through  an  incremental 
area  dA  on  the  wavefront  surface.  While  Eq.  (109)  was  developed  for 
propagating  plane  waves,  it  is  possible  to  arrive  at  the  same  result 
for  spherical  waves  from  an  isotropic  emitter  by  making  plane  wave  ap¬ 
proximations  over  small  elements  on  the  wavefront  surface. 

Let  the  volume  element  dy  in  Fig.  (6)  be  bouoled  by  surfaces  r  - 

Sl'  r  “  S1  +  dS»  €  "  €i*  €  “  €i  +  d€,  v  -  v1#  and  v  -  v1  +  dv,  where 
(r.  €,  v)  is  analogous  to  the  spherical  coordinate  system  (r,  6,  <f>  ) 
centered  at  the  emitter.  If  d€  and  dv  are  sufficiently  sjjb.11,  then  the 
spherical  wavefront  elements  S1  and  S1  +  dS  may  be  approximated  by 
plane  wave  elements.  Following  the  development  in  Chapter  IV,  Eq.  (109) 
may  be  written 

d6dV  “  *av*“  ^  (13*0 

where  d€dv  is  the  element  of  energy  flux  of  the  time-averaged  Poyntlng 
vector  P&v  at  the  wavefront  surface  dA. 

The  program  calculates  relative  power  density  (l.e. ,  normalised  to 
the  power  density  in  free  space  with  no  earth  surface  present)  by  taking 
the  ratio  of  d€dv  in  the  atmosphere  to  d€QdVQ  in  vacuum.  Since  the  at¬ 
mosphere  generally  has  a  refractive  index  gradient  which  varies  with  al¬ 
titude,  Eq.  (126)  predicts  that  the  rays  will  have  differing  radii  of 
curvature  depending  upon  the  local  values  of  ray  slope  ani  Vn/n.  This 
variation  in  ray  curvature  causes  the  rays  to  converge  and  diverge 
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Fig.  6.  The  volume  element  dy  bounded  by  rays  and  surfaces  of 
constant  phase  for  a  spherical  wavefront  element  radiated  from  an  iso¬ 
tropic  emitter  at  point  0.  The  face  abed  lies  on  the  surface  r  -  S1 
and  the  face  efgh  lies  on  the  surface  r  “  S..  +  dS.  Faces  abfe  and  iegh 
are  portions  of  the  surfaces  €  -  €.  and  €  +  d£,  respectively.  The 

remaining  two  faces  adhe  and  begf  lie  on  the  surfaces  v  «■>  w.  and  v  *  v 
+  dv,  respectively.  1  1 


throughout  the  atmosphere,  resulting  in  d€dv  /  d€Qdvo  for  a  given  wave- 
front  element  dA.  However  since  the  atmosphere  is  assumed  to  have  azi¬ 
muthal  symmetry  about  the  emitter,  then  7n  *  0  in  the  v  direction,  which 

yields  dv  -  dvQ.  Thus  the  surface  element  dA  contributes  to  the  flux 

A 

ratio  only  in  the  €  direction.  Replacing  dA  with  its  projection  dl 

,  A 

along  €,  the  flux  ratio  then  becomes 


P 

avo 


■  ndl 


(135) 


Bq-  (35)  is  solved  by  the  following  process,  part  of  which  has  been 
adapted  from  a  ray  tracing  program  developed  by  the  MITRE  Corporation.^ 
At  regular  intervals  along  the  earth's  surface,  which  are  selected  by 
the  user,  the  simulation  divides  the  altitude  difference  between  the 
highest  and  lowest  rays  into  100  slightly  overlapping  vertical  segments, 
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each  of  which  is  used  to  calculate  the  relative  power  density.  Fig.  (7) 
shows  one  such  vertical  (i.e.,  radial  from  the  center  of  the  earth  at  a 
distance  x  from  the  transmitter)  segment  whose  length  is  dly.  The 
angular  subtense  of  the  ray  bundle  which  passes  through  dly  in  free 
space  is  d€Q,  while  the  subtense  of  the  ray  bundle  through  dly  in  a  re¬ 
fractive  medium  is  d€.  Stated  differently,  d€0  and  d€  represent  the 
angular  subtense  of  the  element  dly  as  seen  %  the  emitter  at  point  0 
both  in  free  space  and  in  a  refractive  atmosphere,  respectively. 


EfrgTH 


7.  The  angular  subtense  of  a  line  segment  dl  at  the  emitter 
0  as  seen  in  free  space  and  in  a  refractive  atmosphere . v  Fig.  (7a)  shows 
d€  and  d€  above  a  spherical  earth  surface,  and  Fig.  (7b)  gives  an  en¬ 
larged  view  at  segment  dly. 

Replacing  ndl  in  Eq.  (135)  with  the  line  vector  nydly  gives  the  flux 


ratio  as 
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P  *n  dl 
av  v  v 

P  *n  dl 
a  vo  v  v 


P  dl  cos  € 
av  v 

P  dl  cos  € 
avo  v  o 


P&v  008  € 

’  Pavo  008  €0 


(136) 


where  ny  is  the  unit  vector  normal  to  dly  which  defines  the  local  hori¬ 
zontal  reference  axis,  and  £q  and  £  are  the  respective  elevation  angles 
of  the  Poynting  vectors  ?yo  and  Pay  at  dly.  Solving  Eq.  (136)  for  the 
average  relative  power  density  gives 


P  __  j  £  cos  £ 

p  ,  av  d£ _ o 

avr  P  d€  *cos  £  (137) 

avo  0 

The  angles  d£Q  and  £q  in  Eq.  (137 )  are  obtained  from  the  free  space 
geometry  of  Fig.  (8),  where  a  is  the  radius  of  the  earth,  hQ  and  h  are 
the  known  respective  altitudes  above  the  earth’s  surface  of  the  emitter 
and  a  given  ray,  rQ  and  r  are  the  radial  distances  to  the  emitter  and 
ray  (r0  “  a  +  hQ  and  r  -  a  +  h),  x  is  the  distance  along  the  earth’s 
surface,  and  e  -  x/a  is  the  angle  subtended  by  x  at  the  center  of  the 

earth.  A  right  triangle  is  constructed  having  the  sides  b  and  c  such 
that 


and 


b  ■  rQ  tan6 


(138) 


e  _ _ o_ 

CO80  (139) 

where  side  b  forms  the  local  horizontal  reference  axis  at  the  emitter. 

The  distance  d  between  the  ray  and  emitter  is  given  by  the  law  of  co¬ 
sines  as 


Similarly 


d  ■  [jro2  +  r2  -  2rrocos0]1/2 


f  -  [b2  +  d2  -  2bd  cos  |€ol|  31/2  -|c  -  rl 


(I**)) 

(141) 
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setting  h  =  h,T  and  h  =  h^,  where  h.j  and  ru  are  the  known  upper  and  lower 


height  limits  of  the  vertical  line  segment  dl  in  Fig.  (?).  Thus 


d?o  =  'Uo 


'Lo 


(143) 


ins  angle  ?  between  ?  and  n  is  found  by  reversing  the  posi- 
o  avo  v  ~ 

tions  of  the  emitter  and  ray  in  Fig.  (8)  (i.e.,  by  interchanging  hQ 

with  h  and  r  with  r),  and  replacing  €  .  with  £  .  Since  the  waves  with- 
o  oi  o 

in  the  flux  angular  element  d£o  are  approximately  plane,  the  Poynting 

vector  ?  VQ  in  Fig.  (7)  points  radially  outward  from  the  emitter  and 

through  the  center  of  the  wavefront  surface  element  within  d€Q.  Thus 

setting  h  =  h  in  Eq.  (142)  where  h  is  the  known  height  of  the  center 
c  c 


of  alv,  yields 


€  =  cos 
o 


*[: 


b2  +  d2  -  (c  -  ro)2"l1/2 


2bd 


- 


(144) 


where  b  and  c  are  given  by  2qs.  (138)  and  (139)  with  r  replacing  tq, 
and  d  remaining  unchanged.  ';ith  the  reversal  of  the  emitter  and  ray 
positions  in  Fig.  (8),  the  line  b  now  lies  in  the  direction  of  n^, 
while  line  d  lies  along  ?avo*  Thus  Eq.  (144)  gives  the  angle  between 
the  two  vectors,  both  of  which  are  pointing  to  the  left  in  Fig.  (S). 

The  problem  of  obtaining  the  flux  angle  increment  dc  in  2q.  (137) 
is  different  from  that  of  finding  d€Q,  where  the  elevation  angles  axe 
computed  directly  from  the  ray  altitudes  passing  through  the  upper  and 
lower  limits  of  dl^.  Here  a  simple  geometric  transformation  between 
altitude  and  ray  angle  is  not  possible  since  ray  curvature  is  not  easily 
■odeled  along  each  ray  path.  Instead  a  linear  interpolation  is  made 
between  the  known  initial  elevation  angles  of  the  rays  which  pass  through 
and  on  either  side  of  the  line  segment  dlv. 
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FlS.  9  «  -is ,y  geometry  at  the  vertical  line  segment  dl  within  a 
refractive  atmosphere  v 

Consider  the  set  of  four  rays  shown  in  Fig.  (9)  which  have  initial 
elevation  angles  of  €^f  z^,  and  at  the  emitter.  The  ray 

altitudes  at  the  segment  dlv  axe  ,  h^,  and  h^,  and  the  local  ray 

elevation  angles  are  €^,  ?2,  and  respectively.  The  angles 
and  6^  denote  the  launch  angles  of  rays  passing  through  the  limits  of 
dly  at  heights  h^  and  .  3oth  €T^  and  ST i  are  approximated  by  lin¬ 
early  interpolating  between  the  launch  angles  of  the  rays  which  lie 
immediately  above  and  below  h^T  (rays  1  and  2)  and  those  which  lie  im¬ 
mediately  above  and  below  h,  (rays  3  and  4),  respectively.  The  inter- 

xj 

polation  uses  the  relative  spacing  between  the  rays  and  the  line  seg¬ 
ment  limits  according  to  the  relationships 


”Ji  "li 


-  A? 


,  (hi  :  V 


^  -  h2; 


(145) 


where  (=  0.02  degree)  is  the  equal  angular  separation  between  the 
initial  ray  elevation  angles,  and  dh.j/dh^  and  dh^/dh^  are  the  frac¬ 
tional  differences  in  altitude  between  the  rays  and  the  limits  of  dlv. 
The  flux  angle  for  the  refractive  atmospheric  medium  is  then  computed 


from  Sqs.  (145)  and  (l46)  by 


d~  =  €Ui  "  €Li 


(14?) 


It  should  be  noted  that  Sqs.  (145)  and  (l46)  are  valid  regardless 
of  the  number  of  rays  passing  through  dlv  in  Fig.  (9).  For  instance  if 
rays  2  and  3  were  absent,  the  ratios  for  ?Ui  and  €Li  would  be  = 

€li  "  A'(dhu/dhl4^  and  €Li  =  €li  "  A€(dhi/dhi4^  where  dh*j  =  ^  "  hU’ 
dhL  =  ^1  “  ^L’  and  ddi4  =  ^  The  presence  of  rays  2  and  3  merely 

adds  to  the  accuracy  of  the  interpolations  in  regions  where  the  ray  al¬ 
titudes  are  not  equally  spaced. 

Finally,  the  angle  €  between  F&v  and  ny  is  obtained  by  taking  the 

mean  value  of  the  angles  €^,  and  since  the  Poynting  vector 

is  again  assumed  to  be  pointing  radially  outward  from  the  emitter  along 

the  direction  of  the  rays  and  through  the  center  of  the  flux  element  d€. 

That  is,  the  mean  elevation  an°:le  £  defines  the  elevation  angle  of  P 

av 

at  the  segment  dlv. 

The  elevation  angle  €  of  a  ray  in  polar  coordinates  is  obtained 

ray 

from  Fig.  (2),  where 
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Substituting  P  *  x/a  from  La.  (60b)  and  solving  2q.  (145)  for  €. 


ray 


gives 


'ray 


=  tan 


'  [l  i]  -  ^  [rg  s] 


(149) 


where  h  and  dh/dx  are  solved  in  the  computer  program  by  integrating 
3qs.  (131)  and  (132)  along  the  ray  path.  Using  Eq.  (149)  to  find  the 


elevation  angles  of  the  rays  at  dly,  the  mean  value  of  €  becomes 

c  +  c.  +  c  +  c 


(150) 


which,  like  d€,  improves  in  accuracy  as  the  number  of  rays  passing 
through  dlv  increases. 

Sqs.  (142)  through  (150 )  are  used  in  the  simulation  to  solve  for 
the  relative  average  power  density  in  2a.  (13?)  for  each  of  the  100 
vertical  line  segments  of  length  dlv<  Similarly,  the  program  calcu¬ 
lates  the  magnitude  of  the  relative  field  strength  by  means  of  the  re¬ 
lationship  between  power  density  and  field  strength  shovm  in  £0.  (105). 
Substituting  3qs.  (2),  (3),  (4),  and  (105)  (where  €,  €q  and  ?r  denote 
permittivities)  into  Eq.  (137)  gives 


avr 


n 


d?  .  cos^o 

d€„  cos? 
o 


(151-) 


where  n  is  the  mean  refractive  index  of  air  taken  at  the  center  of  dlv, 
and  3  and  E0  are  the  amplitudes  of  the  electric  field  intensity  vectors 
in  the  refractive  atmosphere  and  in  free  space,  respectively.  Solving 


for  |s/so  | 


gives 


(152) 
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E  r  dg  cosgo  n  1/2 

®r  *  E  ^  n  d€  cos^ 
o  o 

d€  cosS  1/2 

~  r _ 2n 

L  d€  cos<=J 

O 

where  E  is  the  magnitude  of  the  relative  field  strength  and  “v/xT  *  1 
r  v 

in  air. 

If  more  than  one  wavefront  from  a  given  emitter  arrives  at  dly 
because  of  strong  atmospheric  refraction  or  reflection  from  the  earth, 
the  program  computes  the  relative  power  density  and  field  strength  for 
the  resultant  field  and  for  the  field  of  each  component  wavefront  at 
dl  .  Furthermore,  if  the  emitter  radiates  pulsed  rather  than  continu¬ 
ous  waves  (as  in  the  case  of  most  radars)  the  calculations  for  the  re¬ 
sultant  field  include  only  those  component  waves  whose  pulses  overlap 
in  time  with  the  pulses  of  the  direct  (line -of -sight)  propagating  wave. 
Since  the  optical  path  length,  and  hence  the  time  of  propagation,  from 
the  emitter  to  dlv  is  generally  different  for  each  component  wave,  this 
overlap  will  not  extend  over  the  full  pulse  width  of  the  line-of -sight 
wave,  given  that  the  pulse  widths  of  each  wave  are  equal.  Therefore 
all  calculations  for  the  resultant  field  at  dly  denote  maximum  values 
of  Pavr  and  Er  when  all  pulses  overlap  with  the  pulse  of  the  line-of - 
sight  wave. 

Consider  the  case  of  M  waves  originating  from  the  same  emitter 
which  arrive  at  the  vertical  line  segment  dly  at  some  time  t.  Let  the 
emitter  be  linearly  polarized  so  that  the  electric  field  of  each  wave 
will  be  polarized  in  the  same  direction  (to  eliminate  elliptical  and 
circular  polarizations).  Since  Eqs.  (151)  and  (152)  apply  to  those 
field  components  whose  Poyntlng  vectors  are  perpendicular  to  dlv,  the 
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resultant  or  total  field  which  is  propagating  normally  to  dly  may  be 
written  as 

K 

S^e-^  -  2  v.  e“JeiB  (153) 

m-1  m 

where  3  and  6  are  the  amplitude  and  phase  of  each  field  component 
m  m 

propagating  normally  to  dlv. 

For  the  purpose  of  simulation,  all  field  magnitudes  and  phases  are 
represented  by 


a 

m 


and 


(154) 


m 


6+0 

m 


m 


(155) 


where  D  is  the  reflection  coefficient  of  the  earth's  surface.  E  '  is 
/  m  m 

the  field  strength  due  to  atmospheric  refraction,  6^  is  the  wavefront 

chase  at  dl  ,  and  0  is  the  phase  shift  due  to  reflection.  The  wave- 
v  T  m 


front  phase  is  computed  from 


6  ■  kL  -  (k  n)(ct  ) 
m  m  '  o  m' 


(1 56) 


where  the  optical  path  length  L  between  the  emitter  and  dl  is  obtained 

m  v 

by  multiplying  the  speed  of  light  c  in  vacuum  times  the  total  wavefront 

propagation  time  t  as  given  by  Eq,  (133)* 

In  the  simulation  O  ■  1  and  0  *  0  for  each  wavefront  until 

I  m  m 

it  is  reflected  from  the  earth's  surface,  at  which  point  O  and  0 

l  m  m 

are  calculated  as  functions  of  emitter  frequency  and  polarization,  ter¬ 
rain  composition,  surface  roughness,  and  wave  incidence  angle  at  the 
ground.  In  the  case  of  multiple  reflections,  both  terms  are  altered 
accordingly,  with  p ^  equal  to  the  product  of  the  reflection  coeffici¬ 
ent  magnitudes  at  each  reflection  point  and  0m  equal  to  the  sum  of  the 
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reflection  phase  shifts. 

Taking  the  ratio  of  the  total  electric  field,  in  1c.  (153)  to  the 
free  space  electric  field  (without  the  earth's  surface)  gives  the  total 
relative  field 


h,Pe"jeT  r-» 

EU  e  jeTr  =  — - 7T-  =  > 

-r  -n  “j6  / 

.=/  e  o  L i 


O  3  '  e  *^m  + 

I  m  m  _ 


2  e'j6o 
o 


(157) 


m=l 


where  PQ  =  1  and  aQ  =  0  for  the  free  space  field.  Substituting  Eq. 


(152)  for  E  ’/E  in  So.  (157)  gives 

m=l  L 


1/2 


-j(  6  +  <f)  -  b  ) 


e  " '  m  '  m  o 


pl/2 


K 


-KT  ME  I «. [^] 


1/2 


-j(^m  + 


e  wx  m 


(153; 


Since  2_  is  the  magnitude  of  the  total  relative  field  at  dl  ,  the  aver- 
_r  J  v* 

age  relative  power  density  may  be  found  from  Eq.  (151 )  as 


iavr 


(159) 


Eos.  (152)  and  (159),  which  are  merely  extensions  of  Eqs.  (152)  and 

(151),  are  used  in  the  simulation  for  all  relative  field  strength  and 

power  density  calculations.  The  coefficients  O  .  <b  and  6  for  each 

/  m  ~m  m 

reflected  wave  crossing  dl  are  found  in  much  the  same  manner  as  the  wave 

v 

elevation  angle  €  in  Eq.  (150)*  Consider  the  rays  from  a  reflected  wave- 
front  which  pass  through  and  on  either  side  of  dl  as  shown  in  Fig.  (9). 
Since  the  restriction  of  5o.  (115)  in  Chapter  IV  requires  that  adjacent 
raps  be  nearly  parallel  (i.e.,  the  relative  spacing  between  adjacent  rays 


5$ 

must  not  change  appreciably  over  a  wavelength  of  distance),  the  rays 
through  a  small  wavefront  element  may  be  expected  to  intersect  the  , earth 
at  approximately  the  same  angle  of  incidence.  Thus  the  wavefront  re¬ 
flection  coefficients  D  and  .  which  are  dependent  upon  the  angle  of 

l  n  m 

wave  incidence,  may  be  obtained  by  taking  the  mean  of  the  reflection  co¬ 
efficients  P  and  d>  which  are  calculated  at  the  ray  incidence 

i  mp  r  mp 

angles.  For  the  four  rays  in  Fig.  (9),  the  magnitude  of  the  wavefront 

reflection  coefficient  P  becomes 

i  m 

p  e  Pml  P Eft  ,  X  ^2  (l6o) 

4  p«l 

with  similar  expressions  existing  for  4>^  and  the  time  of  propagation  t 


used  in  Sq.  (156)  to  obtain  the  phase  5^.  The  reflection  model  which 


computes  the  values  of  p^  and  is  presented  in  the  following  sec- 


mp 


tion. 


Earth  Reflection  Model 

Consider  a  uniform  plane  wave  which  is  incident  upon  a  plane  bound¬ 
ary  between  two  media  as  shown  in  Fig.  (10).  Let  the  incident  wave  be 
polarized  so  that  its  electric  field  vector  is  either  parallel  to  or 
perpendicular  to  the  plane  of  incidence  (i.e.,  the  plane  containing  the 
wave  propagation  vector  k  and  the  normal  to  the  boundary).  The  ratio  of 
the  reflected  to  the  incident  field  gives  the  reflection  coefficient  y 
which  is  derived  in  any  standard  electromagnetics  text,^'  For  par¬ 
allel  polarization  this  becomes 


*72  cos  fc,  -  7^  cos  ^ 
5“  ‘  Vgi'par  "  ??2  cos  ^  +  7^  cos  ^ 


Sr 

(#) 


(161) 


and  for  perpendicular  polarization 


(162) 


5? 

A  7?2SeC-2  - 

YS“P  ’  Sj  ^  ‘  772secp2  + 

where  Tj  is  the  intrinsic  impedance  of  each  medium,  and  p  denotes  the 
ancle  of  incidence  ani  transmission  for  media  1  and  2  as  shown  in  Fig. 
(10).  For  simple  non-magnetic  media  the  intrinsic  impedance  may  be 


Fig.  10.  Reflection  of  uniform  plane  waves  at  a  plane  boundary  for 
parallel  and  perpendicular  polarizations. 


written  as 


(163) 


where  4,  €,  and  0  are  the  permeability,  permittivity,  and  conductivity  of 
the  medium,  and  are  treated  as  real  quantities.  Note  that  €  is  used  here 
as  a  material  property,  rather  than  an  angular  quantity  as  in  the  preced¬ 
ing  section. 


For  problems  involving  the  reflection  of  a  field  from  the  surface  of 
the  earth,  Eqs.  (l6l)  and  (162)  are  frequently  written  in  a  form  which  is 
dependent  only  upon  the  angle  of  grazing  incidence  0^  where  ■  v/2  - 


5S 

and  the  complex  permittivity  €  of  the  earth  where  €,  which  is  also 

C 

called  the  "a-c  capacitivity”  in  some  texts,  is  the  real  part  of  €  .  To 

C 

derive  these  equations,  which  sire  known  as  Fresnel’s  equations  for  a 
smooth  plane  surface,  it  is  necessary  to  give  the  relationships  between 
the  index  of  refraction,  permittivity,  and  wave  propagation  constant  of 
general  non-magnetic  media  such  as  most  lard  and  sea  surfaces  of  the 
earth.  These  relationships  are  analogous  to  Eqs.  (4)  and  (68)  for  the 
nearly  perfect  dielectric  medium  of  air. 

In  general  both  the  permittivity  and  permeability  of  a  medium  are 
complex  quantities.  However  for  magnetically  lossless  media  such  as 
land  and  sea  bodies,  these  quantities  may  be  given  by 

4C  -  4’  -  *  4*  *  404r  *  40  (l64a) 


€'  -  J€"  «  €c€r  -  j€"  -  €  -  j€" 


(164b) 


where  the  imaginary  term  €"  is  the  dielectric  loss  factor.  Substituting 
Eqs.  (164)  into  Sq.  (4)  gives  the  refractive  index  as 


n 


V 

_E 

c 


-  # 

4o€o 


(165) 


where  ?cr  is  the  relative  complex  permittivity,  or  complex  dielectric 


constant,  and  is  merely  the  complex  form  of  the  relative  permittivity  €r 
or  (real)  dielectric  constant  of  Eq.  (4). 

A  similar  relationship  exists  between  the  complex  dielectric  constant 
and  the  wave  propagation  constant  of  Eq.  (68).  Noting  that  the  general 
-orr.  of  the  wave  constant  is  given  by 


k  -  +  JO>€) 


(166) 


then  taking  the  square  of  Eq.  (70),  solving  for  n,  and  substituting  Eqs. 
(67),  (1 64),  (165),  and  (166)  yields 


€c  ,k  ,2  10  l*o£  - 

C  xz  —  33  [  _  )  *3  .  ■■  I-  —  — — 


"cr  €  vk  ' 
o  o 


g  -  30/C J  _  €*  -  .jg" 


(167) 


where  €"  *  tf/CU  .  Eq.  (167)  is  frequently  expressed  in  mks  units  as 


«cr*f  -  J«°  Xa-e=r'  -  *cr" 
o 


(168) 


where  the  wavelength  X  is  in  meters  and  the  conductivity  0  is  in  mhos/ 
meter. 

Returning  to  the  reflection  coefficient  for  parallel  polarized 
fields  and  letting  the  subscripts  1  and  2  denote  the  air  and  earth  media 
in  Eq.  (l6l),  the  respective  intrinsic  impedances  become 


u  x/2  I*  3 

Valx-  \  -  fe2]  -trs2-] 

aix  1  Vri 


(169) 


since  «  0,  and 


'earth 


laz  +  jO)g2J 


(170) 


Applying  Snell's  law  for  plane  dielectric  media  at  the  earth  surface 
boundary  and  solving  in  terms  of  @2  gives 

hn 

sin  &2  -  (—)  sin  ^ 


vr.ich  may  be  rewritten  as 


(171) 
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n,  2  _  1/2 

cos  £2  =  L  1  -  (~)  Sin  ^  ] 


(172) 


Substituting  2qs.  (4)  and  (165),  for  media  1  and  2  respectively,  into  2q. 
(172'  ~ives 


cos  e9  -  L 1  -  feS-) 


rl  s  .  2  e  n 
- )  sin  £. 


(173) 


Multiplying  the  numerator  and  denominator  of  2o.  (l6l)  by  V1/V22  £ives 
the  reflection  coefficient  for  parallel  polarized  fields  as 

»  °°S  ^  '  (7?l/7?2)2  C0£  %  (m) 

par  cos  e2  +  cn1/V2)2  cos  h 

Noting  that  the  impedance  ratio  at  the  earth  surface  boundary  is 


Vi  f  “o  c2  * Ja) -2I 

^"[Vr!  ^»0  J 


'2  •  jc2/« 
€ 
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*'cr2/rl 


(175) 


then  substituting  2qs.  (173)  and  (174)  into  Eq.  (175)  yields 

^  "ll/2 


’4«)  -  sin2  6,1  -  £*)  coS  (L 

L  rl _  ~rl  x 


(f2^)  -  sin^  g 
-  rl  . 


'2  e  _ 

+  (t3^)  cos  e. 
rl  x 


(176) 


rollowlng  the  same  type  of  procedure  for  Eq.  (162),  the  reflection  coef¬ 
ficient  for  perpendicular  polarized  fields  becomes 

T  9  1^/2 


cos  fc. 


■>  ja£  -  -2 1' 

h  * 


11/ 


(177) 
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Table  (l),  which  is  taken  from  Xerr,^  gives  typical  values  of  the 

real  €  '  and  imaginary  € _ "  terms  of  the  complex  dielectric  constant 

cr  cr 

for  various  water  and  soil  types  at  several  different  frequencies. 

While  these  values  are  not  intended  to  be  a  complete  set  of  earth  para¬ 
meters,  they  do  represent  values  obtained  from  a  variety  of  sources 
which  serve  as  a  useful  guide  for  the  complex  dielectric  constant  of 
earth  surfaces. 


TA3LE  1 


APPROXIMATE  ELECTROMAGNETIC- PROPERTIES 
OF  SOIL  AND  WATER1 


Medium 

X 

G 

ohm/m 

c  * 
"cr 

€  ” 
cr 

Sea  water 

3  m 

4.3 

30 

77  4 

20°  -  25°  C 

20  cm 

4.3 

80 

52 

10  cm 

6.5 

69 

39 

28°  C 

3.2  cm 

16 

65 

30.7 

Distilled  water,  23  C 

3.2  cm 

12 

67 

23 

Fresh  water  lakes 

1  m 

0.001 

80 

0.06 

1  m 

0.01 

80 

0.60 

Very  dry  sandy  loam 

9  cm 

0.03 

2 

1.62 

Very  wet  sandy  loam 

9  cm 

0.6 

24 

32.4 

Very  dry  ground 

1  m 

0.0001 

4 

0.006 

Moist  ground 

1  m 

0.01 

30 

0.6 

Arizona  soil 

3.2  cm 

0.10 

3.2 

0.19 

Austin,  Tex.  soil, 
very  dry 

3.2  cm 

0.0074 

2.8 

0.014 

As  Table  (l)  indicates,  the  properties  of  sea  water  (€  '  and  0 )  are 

cr 

independent  of  wavelength  for  wavelengths  greater  than  20  cm.  Since  €  " 

cr 

*  c/CJ  €q  «  6 0  GX  from  Eq.  (l68),  then  is  found  to  vary  only  with 

wavelength  and  not  sea  conductivity.  However  at  shorter  wavelengths  G 

is  dependent  on  wavelength,  which  affects  the  value  of  € _ "  as  shown  in 

cr 

?ig.  (11).  Furthermore,  €  '  and  G  are  highly  dependent  on  temperature 

cr 

^cr'  ^ecreasing  and  o  increasing  with  higher  temperatures.  Further 


evidence  of  the  temperature  dependence  of  and  €  "  far  sea  and 

33  34 

fresh  water  is  given  "by  Burrows  and  Attwood  and  by  Saxton. 


Fig.  11.  VJavelength  dependence  of  €  '  and  €  "  far  sea  water  at 

O  ^ 

17  C,  from  burrows  and  Attwood.  The  dotted  line  for  £  "  is  far  a 

cr 

constant  0*  3.51  aho/m,  while  the  solid  €  '*  curve  results  from  the 

cr 

dependence  of  C  on  wavelength  for  wavelengths  shorter  than  20  cm. 


For  land  surfaces  both  €  *  and  (J  are  much  lower  than  for  water, 

with  the  smaller  values  associated  with  dry,  rocky,  or  sandy  soil,  and 

the  higher  values  occurring  with  moist  and  rich  soil.  As  Kerr  points _ 

out  in  his  discussion  of  the  values  in  Table  (l),  the  wide  range  in  €  ' 

cr 

and  0  produces  a  considerable  variation  in  €  which  in  turn  largely  ef- 

cr 

iects  the  reflection  coefficient,  especially  for  parallel  polarization. 
*n  tie  other  hand,  the  reflection  coefficient  for  either  polarization  is 


not  appreciably  affected  even  if  medium  1  is  assumed  to  be  free  space 
* 'ri  *  l)i  since  the  mean  relative  permittivity  of  air  at  sea  level  is 


(178) 


€r(aix)  *  €rl  "  “l2  *  U-00031)2  «  1 

and  Jc^j  >  1  for  extremely  dry  soil  and  |  »  1  for  all  other  soil 

and  water  surfaces  in  Table  (l). 

Frequently  the  coefficient  of  reflection  is  described  in  terms  of 

the  grazing  incidence  angle  shown  in  Fig.  (10 )  as  ■  tt/2  and  the 

polarization  of  the  electric  field  with  respect  to  the  earth's  surface 

where  the  terms  "vertical"  and  "horizontal"  polarization  axe  synonymous 

with  parallel  and  perpendicular  polarization.  Taking  note  that  €  ,  «  1. 

rl 

Eqs.  (l?6)  and  (177)  may  then  be  written  as 


^or2  ~  cos2  °1^2  *  ecr2  sl°  S 
t€cr2  *  0062  “l^2  +  ser2  sln  °i 


(179) 


and 


0  .a. -A  -  - -  5.  cVz ' 

n  •  Jl  «  , 

sin  0^  +  [€^2  -  cos"  0^ 


-  cos" 

T 


f/2 

J7z 


(180) 


where  the  subscripts  v  and  h  denote  vertical  (parallel)  and  horizontal 
(perpendicular)  polarizations,  and  p  and  <t>  are  the  magnitude  and  phase 
of  the  reflection  coefficients. 

Figs.  (12)  through  (17),  which  are  taken  from  Long^-5  and  Povejsil, 

36 

Raven,  and  Waterman,  show  the  amplitude  and  phase  of  the  Fresnel  re¬ 
flection  coefficient  of  Eqs.  (179)  and  (l80)  for  both  vertical  and  hori¬ 
zontal  polarizations.  Figs.  (12)  through  (15)  are  for  a  smooth  sea  sur¬ 
face  at  10°  C  with  €  ’  and  C  having  values  comparable  to  those  of  Table 

(l).  The  most  pronounced  feature  is  the  insensitivity  of  reflected  hori¬ 
zontally  polarized  fields  to  grazing  angle  as  compared  to  vertically 
polarized  fields.  For  vertical  polarizations,  the  reflection  coefficient 
magnitude  is  at  a  minimum  when  -  tt/2.  The  corresponding  grazing 


14,  Phase  of  the  reflection  coefficient  as 
grasing  ancle  for  a  snooth  sea  at  10°  C  ^ 


15.  -expanded  plot  of  ?i~. 
and  10  decrees 


14  for  incident  grazing  angles 
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?ig.  1 6.  Magnitude  of  the  reflection  coefficient  as  a  function 
of  incident  grazing  angle  for  smooth  average  land  where  €'  ■»  10  and 
C  *=  .0016  nho/n  35,  3°  cr 
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1?»  Phase  of  the  reflection  coefficient  as  a  function  of 
incident  grazing  angle  for  smooth  average  land  where  €'  ■  10  and  o 
■  .COlo  mho/m  35,  35 


6? 


ancle  a  is  known  as  Brewster's  ancle  which  represents  the  angle  of  in¬ 
cidence  for  which  transmission  into  the  sea  is  maximized.  Since  the 
impedance  ratio  of  2q.  (175)  is  V^Vz  =  "V^cr2^1  ~  "N/  ?cr2  and 
since  j  -c„2  J  >;>  i  for  sea  water,  the  impedance  mismatch  at  the  air 
and  sea  boundary  will  always  result  in  wave  reflections  even  at  the 
Brewster  angle. 

As  previously  mentioned,  the  values  of  €  '  and  c  are  much  smal- 

C2T 

ler  for  land  surfaces  than  for  sea  or  fresh  water  surfaces.  Land  is 
generally  a  better  dielectric  material  than  water  and  yields  a  consid¬ 


erably  improved  Impedance  match  at  the  surface.  This  is  evident  in 
Pigs.  (l6)  and  (17 )  where  the  values  =  '  =10  and  c  =  0.0016  r.ho/n 

are  used.  Pig.  (lo)  shows  that  p  is  iiearly  zero  for  vertically  polar¬ 
ized  fields  incident  on  smooth  average  land  at  the  Brewster  angle. 

This  contrasts  with  values  of  approximately  p  =  0.1  to  p  -  0.4  at 
the  Brewster  angle  for  smooth  sea  surfaces.  Furthermore,  the  smaller 
value  of  J'^J  results  in  a  general  decrease  in  p  for  both  hori¬ 
zontal  and  vertical  polarizations  and  a  nearly  instantaneous  phase 
change  at  Brewster's  angle  for  vertically  polarized  fields. 

Since  the  earth’s  electrical  properties  vary  so  considerably,  sev¬ 
eral  surfaces  are  modeled  in  the  simulation  by  using  representative 
values  of  ?  '  and  o  as  shown  in  Table  (2).  The  sea  models  are  divided 

into  three  wavelength  regions  whereas  the  land  models  are  relatively 
independent  of  waveleagth  but  dependent  upon  the  moisture  content  of 
*.".e  ground.  All  values  are  averages  obtained  from  a  number  of  sources 
end  are  therefore  only  approximate.^-'  ^ 

'-n  reality,  the  Fresnel  coefficients  are  of  limited  use  in  describ- 
-ft.  reflected  fields  since  few  earth  surfaces  may  be  regarded  as  being 
ei.r.cr  plane  or  smooth.  Consequently,  a  modification  to  £qs.  (179)  and 
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TABLE  2 


ELECTROMAGNETIC  PROPERTIES 
OF  SOIL  AND  WATER  MODELS 


Medium 

X 

c 

mho/m 

c  » 

"CT 

3ea  water 

10 

m 

_ 

20 

cm 

4.3 

80.0 

20 

cm 

- 

6 

cm 

6.5 

69.0 

6 

cm 

- 

1 

cm 

16.0 

65.0 

Met  ground 

10 

m 

- 

1 

cm 

.01 

30.0 

Average  ground 

10 

m 

- 

1 

cm 

.0016 

10.0 

Dry  ground 

10 

m 

*■ 

1 

cm 

.0001 

4.0 

(ISO)  is  required  for  most  applications.  Kerr,^  Beckmann  and  Spizzi- 
17  18 

chino,  and  Barton  present  such  a  theory  which  is  summarized  in  the 
following  discussion. 


Reflection  from  any  generalized  surface  is  usually  described  in 
terms  of  specular  and  diffuse  scattering,  vihere  "specular"  scattering 
means  that  the  angle  of  incidence  equals  the  angle  of  reflection  and 
"diffuse"  scattering  implies  reflection  in  all  other  directions.  How¬ 
ever  diffuse  scattering  is  a  highly  complicated  function  of  incidence 
angle,,  surface  electrical  properties,  surface  roughness,  wavelength,  and 
geometry.  As  such  it  does  not  lend  itself  to  a  simple  expression  which 
is  valid  for  all  angles  and  wavelengths,  and  therefore  lies  beyond  the 
scope  of  this  thesis.  Specular  reflection,  on  the  other  hand,  is  the 
dominant  type  in  most  applications  and  is  also  readily  solvable. 

The  coefficient  of  specular  reflection  for  the  earth's  surface  is 
generally  given  by 

y  =pe-J$*  yv^d:;  =  pVjhDRe'j^,h  (181) 

";,9ro  Y^y.  “  pv  ^^,h  is  the  Fresnel  coefficient  for  vertical  and 
horizontal  polarizations,  D  is  a  divergence  factor  which  describes  the 
reduction  in  reflection  caused  by  the  earth's  curvature,  and  F:  is  a 
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factor  of  surface  roughness-  D  has  values  which  range  between  1  and  0 
which,  depending  on  the  geometry  between  the  emitter,  earth,  and  re¬ 
ceiver,  correspond  to  the  cases  in  which  the  earth  may  be  regarded  as 
flat  or  highly  curved.  The  roughness  factor  2  likewise  ranges  from  1 
to  0  depending  on  whether  the  surface  is  perfectly  smooth  or  extremely 
rough. 

The  geometric  spreading  of  wavefronts  due  to  reflection  from  a 
spherical  earth  surface  is  illustrated  in  Fig.  (l8),  This  wave  diver¬ 
gence  is  already  accounted  for  in  the  simulation  since  all  ray  trajec¬ 
tories  are  calculated  in  polar  coordinates  and  are  reflected  from  a 
spherical  earth.  Thus  the  total  specular  reflection  coefficient  used 
in  the  simulation  becomes 

Y=pe_^pv  hRe“M,h  (182) 


Fig.  18.  Geometry  of  the  divergence  factor  D,  from  Long. ^  The 
dotted  perimeter  5  represents  the  cone  of  rays  having  an  initial  angu¬ 
lar  separation  of  £1  which  are  reflected  from  a  tangent  plane  at  the 
center  of  the  cone.  The  solid  perimeter  3'  shows  the  additional  ray 
divergence  due  to  reflection  from  the  curved  earth  surface. 
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Unlike  divergence,  surface  smoothness  is  considerably  more  diffi¬ 
cult  to  quantify.  In  general  a  smooth  surface  is  defined  as  one  which 

17  35 

satisfies  the  condition  established  by  Lord  Rayleigh  given  by  ' 


i  X 

n  8  sina 


(183) 


where  6h  is  the  height  of  surface  irregularities,  X  is  the  wavelength, 
and  a  is  the  incident  grazing  angle  between  the  ray  and  a  plane  surface 
representing  the  average  of  the  surface  irregularities.  Some  authors 
replace  the  factor  8  with  16  in  2q.  (183)  due  to  the  difficulty  of 
clearly  defining  smooth  versus  rough  surfaces,  nonetheless,  the  .tay- 
leigh  criterion  states  that  the  variations  in  height  must  become  smal¬ 
ler  as  wavelength  decreases  and  grazing  angle  increases ,  for  a  surface  .. 
to  be  considered  smooth.  For  example,  if  X  *  1  m  then  6h  must  be  less 
than  approximately  7.1  m  at  a  grazing  angle  of  1  degree,  and  less  than 
1.4  m  at  a  grazing  angle  of  5  degrees.  If  X  *  1  cm  then  the  maximum 
values  of  6h  become  7.1  cm  and  1.4  cm  for  grazing  angles  of  1  and  5  de¬ 
grees  respectively.  Thus  a  calm  sea  might  be  expected  to  appear  as 
more  of  a  smooth  surface  than  most  land  for  wavelengths  ranging  from 
the  radio  spectrum  to  the  microwave  region. 

In  practice  it  is  found  that  the  earth' s  surface  may  usually  be 
described  in  terms  of  a  Gaussian  distribution  of  surface  heights  about.- 
s.  local  mean  height.^*  ^  Ament‘S  and  jeckmann  and  Splzzichino^ 
have  calcul?ted  the  mean  square  value  of  the  specular  roughness  factor 
■'  for  a  Gaussian  surface  (neglecting  shadowing  and  sharp  edge  diffract¬ 


ion  effects)  as 


<  IT  >  -  e 


-U+Y 


(184) 


_  4n£h.slm 
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(185) 


where  <  >  denotes  the  mean  value  of  a  function, Ah  is  the  standard  devi¬ 
ation  of  the  Gaussian  distribution  of  heights,  a  is  the  incident  grazing 
angle,  and  X  is  the  wavelength  of  the  incident  field.  Fig.  (19)  shows 
the  relationship  of  <  >  to  A$,  with  2q.  (184)  plotted  as  a  solid  line 

and  experimental  data  from  sea  and  land  surfaces  (both  plane  and  hilly 
conditions)  represented  by  crosses.  Spizzichino  states  that  the  corre- 
lation  of  <  >  with  A<b  is  good  considering  that  the  reflection  data 

were  often  given  with  little  precision  and  the  values  used  for  Ah  were 
merely  estimated. 

Frequently  land  and  sea  surface  profiles  are  reported  only  in  terms 
of  the  largest  observed  height  variations,  rather  than  as  a  distribution 
of  height  measurements.  The  data  in  Tables  (3)  and  (4),^’  ^  while  rep¬ 
resenting  surface  descriptions  that  can  be  regarded  as  loosely  quantita¬ 
tive  at  best,  often  provide  the  only  basis  available  from  which  Ah  may  be 
obtained.  For  example,  the  Douglas  sea  scale  listed  in  Table  (3)  defines 
wave  height  as  the  average  of  the  peak-to-txough  heights  of  the  one-third 
highest  waves  in  a  given  observation.  Even  if  the  one-third  highest  waves 
are  predominant,  there  is  rarely  any  statistical  measure  established  for 
smaller  waves  in  most  observations,  thereby  rendering  an  accurate  meas¬ 
urement  of  Ah  impossible. 

Similar  difficulty  exists  in  establishing  Ah  for  land  surfaces  where 
the  irregularities  easily  become  an  order  of  magnitude  greater  than  those 
of  the  roughest  seas.  Table  (4)  gives  the  average  of  the  peak-to-mean 
height  variations,  or  deviations  from  the  mean,  of  the  largest  scale 
surface  features  for  several  types  of  terrain.  Land  roughness  is  often 
described  in  terms  of  the  surface  height  deviation  from  the  local  mean 


^ig.  19.  Kean  square  <  >  of  the  surface  roughness  factor  B 

phase  deviation^,  from  Beckmann  and  Spizzichino^? , 


TA3LE  3 

SEA  STATE  AND  WAVE  HEIGHT35 


Douglas  Sea  State 

Description 

Peak-to-Trough 
Have  Height 
(feet) 

Have  Height 
Deviation  Ah 
(feet)  m 

1 

Smooth 

0-1 

0  -  .5 

2 

Slight 

1-3 

.5  -  1.5 

3 

Moderate 

3-5 

1.5  -  2.5 

4 

Rough 

5-3 

2.5  -  4.0 

5 

Very  rough 

8-12 

4.0  -  6.0 

6 

High 

12  -  20 

6.0  -  10.0 

7 

Very  high 

20-40 

10.0  -  20.0 

8 

Precipitous 

40  -  60 

20.0  -  30.0 

TA3LE  4 

LAND  SURFACE  DEVIATIONS39 

Description 

Surface  Height 
Deviation  Ah 
(feet)  m 

Very  smooth  plains 

0-20 

Smooth  plains 

20  -  70 

Slightly  rolling  plains 

70  -  130 

Rolling  plains 

130  -  260 

Hills 

260  -  500 

Mountains 

500  -  1000 

Rugged  mountains 

1000  -  3000 

Extremely  rugged  mountains 

3000  and  above 

surface  altitude  (Ah^  *  h  -  <  h  >)  because  of  the  asymmetrical  nature  of 
terrain  irregularities. 

Given  that  most  earth  surface  profiles  are  approximately  Gaussian, 

the  data  of  Tables  (3)  and  (4)  may  be  Interpreted  as  upper  values  of 

their  respective  height  distributions.  Bullington^  and  other  auth- 
17  37 

ors  ’  suggest  letting  the  maximum  height  deviation  Ah,,  of  a  surface 
profile  represent  the  amplitude  of  the  surface  irregularity  which  is  ex¬ 
ceeded  less  than  1  per  cent  of  the  time  over  the  path  between  the  emitter 
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and  receiver.  Using  this  criterion,  the  maximum  surface  deviations  in 
Tables  (3)  and  (4)  Kill  include  93  per  cent  of  all  height  variations  (49 
per  cent  above  and  below  the  mean  surface  height). 

It  may  be  shown  from  any  set  of  statistical  tables  that  93  per  cent 
oi  all  possible  random  variables  having  a  standard  normal  or  Gaussian 
distribution  will  be  found  within  2.33  standard  deviations  of  the  mean. 
Thus  letting  Ahjj  be  the  maximum  surface  deviation  about  the  mean  such 
that  Ahjj  is  exceeded  less  than  1  per  cent  of  the  time,  then 
Ahjj  «  hM  -  <  h  >  «=  2.33Ah 

where  h^  is  the  maximum  surface  height  in  any  observation.  Values 
Ah^  and  Ah  used  in  the  simulation  are  given  in  Tables  (5)  and  (6), 

Ahjj  is  chosen  as  the  maximum  height  deviation  listed  in  Tables  (3) 

(4)  for  each  surface  description. 


TABLE  5 

STANDARD  DEVIATION  OF  HEIGHTS  FOR  SEA  SURFACES 


Douglas  Sea  State 

Description 

Ah.. 

(feet) 

Ah 

(feet) 

1 

Smooth 

.5 

.2 

2 

Slight 

1.5 

.6 

3 

Moderate 

2.5 

1.1 

4 

Rough 

4.0 

1.7 

5 

Very  rough 

6.0 

2.6 

6 

High 

10.0 

4.3 

7 

Very  high 

20.0 

8.6 

8 

Precipitous 

30.0 

12.9 

Combining  Eqs,  (179)  through  (l36)  and  Tables  (2),  (5),  and  (6),  the 
total  specular  reflection  coefficient  becomes 

Y  p^hRe-H,h  -  py  j y<  R2  >  e'H,h 


(186) 

of 

where 

and 


(187) 
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TA3LE  6 


STANDARD  DEVIATION  OF  HEIGHTS  FOR  LAND  SURFACES 


description 

AJx. 

( feet ) 

Ah 

(feet) 

Very  smooth  plains 

20 

9 

Smooth  plains 

70 

30 

Slightly  rolling  plains 

130 

56 

Rolling  plains 

260 

112 

Hills 

500 

214 

Mountains 

1000 

429 

Rugged  mountains 

3000 

1288 

Extremely  rugged  mountains 

50  00 

2146 

where  E  =  R,  zs“\f<  >.  Because  of  the  assumed  Gaussian  nature  of  the 
rms 

surface  profile,  the  Fresnel  reflection  phase  ^  h  is  the  only  explicit 
term  in  Eq.  (187)  which  contributes  to  a  phase  shift  in  the  reflected 
field.  Hhile  a  given  surface  deviation  will  alternate  the  reflected  field 
by  an  equal  amount  regardless  of  whether  it  lies  above  or  below  the  mean 
surface  height,  the  difference  in  path  length,  and  hence  the  phase  shift, 
between  fields  reflected  from  the  mean  and  the  surface  variation  will 
differ-  in  sign  according  to  whether  the  height  deviation  lies  above  or 
below  the  mean.  Thus  the  phase  changes  due  to  equally  sized  surface  de¬ 
viations  above  and  below  the  mean  surface  height  will  be  equal  but  oppo¬ 
site  in  sign.  Assuming  a  Gaussian  surface  therefore  results  in  random 
phase  shifts  having  a  zero  mean  value  which  do  not  contribute  to  the  phase 
of  the  total  reflection  coefficient  of  Eq.  (187). 


CHAPTER  VI 


RESULTS 

Introduction 

In  an  elimination  of  over-the-horizon  radio  propagation,  Pappert 
lU 

and  Goodhart  recently  presented  an  excellent  comparison  of  theoreti¬ 
cal  and  experimental  results  for  long  range  tropospheric  propagation 
due  to  ducting  conditions  off  the  San  Diego,  California,  coast.  A 
ground  based  duct  was  analyzed  by  means  of  traveguide  mode  theory  and 
field  strength  measurements  at  65,  170,  520,  and  3300  MHz  using  vari¬ 
able  emitter  and  receiver  heights,  while  an  elevated,  o-r  earth  detached, 
duct  iras  studied  at  the  single  frequency  of  3087.7  KHz. 

The  data  presented  by  Pzppert  and  Goodhart  were  used  to  point  out 
the  relative  merits  of  using  waveguide  concepts  to  describe  anomalous 
propagation.  A  similar  comparison  Is  made  in  this  thesis  to  evaluate 
the  ability  of  geometric  optics  to  predict  the  field  In  a  layered  at¬ 
mospheric  structure  such  as  a  duct.  However,  validation  of  the  thesis 
geometric  optics  model  is  restricted  to  the  case  studies  of  the  ground 
based  duct  since  Pappert  and  Goodhart  presented  results  for  the  elevated 
duct  that  were  limited  to  a  single  set  of  transmitter  and  receiver  heights 
and  only  one  frequency. 

The  presence  of  tropospheric  ducts  off  the  California  coast  is  a 
!fs~  fcnown  phenomenon  that  is  largely  caused  by  strong  temperature  in¬ 
versions  ranging  in  height  from  near  sea  level  to  4000  feet.^'  ^3  The 
air  below  these  inversions  is  usually  moist  and  well  mixed  while  the  air 
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above  Is  warmer  and  much  drier.  A  rapid  decrease  in  the  refractive  in¬ 
dex,  as  indicated  by  Eq.  (15),  occurs  in  a  thin  layer  between  these  two 
contrasting  types  of  air,  and  results  in  the  kind  of  ducting  shown  in 
the  Guadalupe  Island  case  studies. 

The  experimental  data  considered  in  this  thesis  were  obtained  from 
field  strength  measurements  taken  by  the  U.  S.  Naval  Electronics  Labora¬ 
tory  Center  across  a  280  nautical  mile  over-water  path  between  San  Diego 
and  Guadalupe  Island.  A  receiver  located  at  heights  of  100  and  500  ft 
above  mean  sea  level  at  San  Diego  recorded  the  signal  of  a  horizontally 
polarized  airborne  transmitter  which  was  flown  in  and  out  of  a  ground 
based  duct  lying  between  these  two  locations. 

The  measured  data,  which  are  shown  in  Fig.  (20)  for  the  100  ft  re¬ 
ceiver,  are  plotted  as  the  received  field  strength  normalized  to  the  free 
space  field  (in  decibels)  versus  transmitter  height.  These  data,  com¬ 
monly  referred  to  as  height  gain  curves,  are  given  at  20  naut  mi  inter¬ 
vals  (measured  along  the  surface  of  the  earth)  between  the  emitter  and 
receiver.  Also  shown  is  the  approximate  flight  path  of  the  transmitter 
and  the  meteorological  profiles  of  atmospheric  refractlvity.  The  dis¬ 
tance  from  the  San  Diego  receiver  to  the  geometric  horizon  is  given  with 
each  set  of  height  gain  curves  to  illustrate  that  the  measured  field  re¬ 
sults  from  over-the-horizon  propagation  due  to  tropospheric  ducting. 

Although  the  data  in  Fig.  (20)  were  obtained  in  19^8 ,  they  represent 
some  of  the  best  case  studies  of  field  measurements  and  supporting  mete¬ 
orological  data  that  are  available  today. ^  The  meteorological  profiles 
shown  along  the  San  Diego  to  Guadalupe  Island  path  are  recordings  of  the 
atmospheric  refractive  index  structure  given  in  3-units  (Chapter  II). 

*hese  profiles  are  easily  converted  to  the  more  conventional  refractivity 


Jefractlvlty  profiles  and  height  gain  curves  for  the  Guadalupe  Island  ground  based  duct 


(133; 


79 

or  N-units  by  combining  Eqs.  (7)  and  (8)  to  yield 


An  examination  of  Fig.  (20)  indicates  that  the  layer  structure  varied 
temporally  or  spatially  both.  I.'one  the  less,  a  strong  gradient  of 
-22  B-units/lOOO  ft  (-22  3-units/kft)  or  -232  N-units/kft  exists  from 
approximately  600  to  1000  ft  in  altitude,  which  exceeds  the  -47. 85  i«- 
units/kft  (-157  N-units/km)  criterion  for  ducting  given  by  Eq.  (10). 

An  a.verage  3-profile,  shown  in  the  lower  left  hand  corner  of  Fig.  (20) 
and  listed  in  Table  (7),  is  used  in  the  thesis  simulation  to  represent 
the  refract ivity  profile  along  the  entire  San  Diego  to  Guadalupe  Island 
path. 

A  tril inear  approximation  to  the  Guadalupe  Island  refractivity  pro¬ 
file  was  similarly  used  in  the  waveguide  computations  of  Pappert  and 
Goodhart.  The  trilinear  model,  which  was  originally  given  in  modified 
index  or  M-units ,  is  converted  to  N-units  by  combining  Eqs.  (?)  and  (9) 

such  that 


and  is  listed  in  Table  (8).  3oth  profiles  are  used  in  the  simulation 
with  the  assumption  that  the  atmospheric  refractivity  structure  is  rea¬ 
sonably  stationary  and  homogeneous  in  the  horizontal  direction. 

Comparisons  of  the  two  profile  models  and  their  N-gradients  are 
given  in  Figs.  (21)  and  (22).  3oth  models  have  extremely  large  negative 
gradients  from  600  to  1000  ft  in  height  and  positive  gradients  between 
1000  and  1300  ft,  which  correspond  respectively  to  superrefractive  and 
subrefractive  regions.  It  is  this  superrefractive  layer  and  the  sea  sur¬ 
face  which  act  as  waveguide  walls  tc  trap  horizontally  traveling  waves, 
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TA3LE  7 


ATMOSPHERIC  REF.JICTIVITY  PROFILE  FOR  THE 
GUADALUPE  ISLAIO  DUCT 


Height  above 
mean  sea  level 


an 


Atmospheric  Refractivity 
B-units  ii-units 


0 

341.0 

341.0 

100 

341.0 

338.8 

200 

341.0 

337.6 

400 

341.0 

336.2 

600 

342.0 

334.8 

700 

320.0 

311.6 

800 

305.0 

295.4 

900 

296.0 

285.2 

1000 

287.0 

275.0 

1100 

294.0 

280.8 

1200 

298.0 

283.6 

1300 

300.0 

234.4 

1450 

301.0 

283.6 

1600 

300.0 

280.8 

1700 

299.0 

278.6 

1300 

298.0 

276.4 

1900 

296.0 

273.3 

2000 

293-0 

269.1 

2100 

294.0 

268.9 

2200 

296.0 

269.7 

2300 

297.0 

269.5 

2500 

297.0 

267.1 

TABLE  8 

ATMOSPHERIC  REFRACTIVITY  .PROFILE  FOR  THE 
TRILINEAR  DUCT 


Height  above 
mean  sea  level 


Ml 


Atmosp:  eric  Refractivity 
H-units  H-units 


.'hile  the  tril  inear 


3] 

thus  enabling  wave  propagation  over  the  horizon, 
model  is  a  good  approximation  of  the  Cuadalupe  Island  duct,  it  under¬ 
estimates  the  strong  N-gradient  existing  from  30C  to  1000  ft.  Fortu¬ 
nately  this  difference  is  not  crucial  since  both  models,  whose  maximum 
gradients  at  this  layer  are  -232  I.-units/kft  (Guadalupe  Island)  and 
-153  N-units/kft  (trilinear),  satisfy  the  -47. 85  H-units/kft  gradient 
required  for  ducting.  A  second  layer,  extending  from  1300  to  2000  ft 
in  the  Guadalupe  Island  profile,  is  shown  to  be  highly  refractive  (-42 
N-units/kft)  although  there  is  no  evidence  in  either  the  height  gain 
curves  of  Fig.  (20)  or  in  the  simulation  results  that  this  layer  causes 
any  significant  trapping. 

The  theoretical  calculations  of  Pappert  and  Goodhart  are  based  upon 

the  modal  wave  solution  for  a  planar  waveguide  as  developed  by  rud- 
14  i« 

den.~  '  ~  In  this  solution  the  transverse  electric  (IS)  wave  propagat¬ 

ing  in  a  trilinear  medium,  such  as  the  one  in  Fig.  (21),  is  obtained 
from  the  plane  wave  reflection  coefficients  at  the  boundaries  of  the 
medium.  The  reflection  coefficients  are  given  as  functions  of  the  eigen¬ 
values  of  the  incident  grazing  angle  for  each  propagating  mode  at  the 
medium  boundaries.  2ach  reflection  coefficient  is  then  expressed  in 
terms  of  modified  rlankel  functions  of  order  one  third  and  their  deriva¬ 
tives,  Entering  into  these  coefficients  are  the  Fresnel  reflection  co¬ 
efficient  and  surface  roughness  factor  given  by  3qs.  (170)  through  (l35), 
In  all  of  their  calculations  Pappert  and  Goodhart  use  a  surface  height 
standard  deviation  (Ah)  of  1  ft,  which  corresponds  to  a  sea  state  of  3 
in  Gable  (5).  Furthermore,  the  waveguide  calculations  require  the  solu¬ 
tion  of  from  one  to  a  hundred  modes  depending  on  the  height  of  the  duct¬ 
ing  layer,  the  frequency,  and  the  emitter  location  within  the  ducting 


Fig.  21.  Conqparison  of  atmospheric  refractivity  profiles 
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Fig.  22.  Comparison  of  atmospheric  refractivity  gradient  profiles 
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Simulation  Results 

Results  obtained  from  the  geometric  optics  computer  simulation  de¬ 
scribed  in  Chapter  V  are  presented  in  this  section  far  both  the  Guadalupe 
Island  and  trilinear  ground  based  duct  refractivity  profiles.  A  moderate 
sea  surface  (sea  state  3)  is  used  to  approximate  the  surface  model  used 
by  Pappert  and  Goodhart.  A  perfectly  smooth  sea  could  have  been  assumed, 
however,  since  the  surface  roughness  factor  of  £q.  (lS4)  is  very  nearly 
unity  (0.959)  at  the  highest  frequency  (3300  i£iz)  and  largest  grazing 
angle  (0.359  deg)  encountered  in  this  study.  The  emitter  characteristics 
are  the  same  as  far  the  Guadalupe  Island  measurements,  namely  a  horizon¬ 
tally  polarized  radiating  pattern  at  65,  170,  520,  and  3300  MHz. 

Since  the  atmosphere  is  considered  to  be  a  linear  and  isotropic  med- 

31 

ium,  the  theorem  of  reciprocity  may  be  applied.  Reciprocity  states 
that  for  a  linear  and  isotropic  medium  the  transmitter  and  receiver  may 
be  interchanged  without  affecting  the  response  of  either.  Thus  the 
transmitter  and  receiver  positions  are  reversed  in  the  simulation  with 
the  emitter  located  at  100  and  500  ft  altitudes  and  the  receiver  moving 
vertically  into  and  out  of  the  duct  along  the  flight  profile  of  Fig.  (20). 

3efore  examining  the  simulation  results  for  propagation  within  a 
ground  based  duct,  consider  the  case  of  an  isotropic  emitter  in  free 
space  (vacuum)  which  is  at  a  height  hQ  of  100  ft  above  mean  sea  level.- 
Fig.  (2 3)  shows  the  refractivity  and  refractivity  gradient  profiles  of 
free  space,  which  are  of  course  zero  since  the  index  of  refraction  n  in 
3q.  (?)  is  unity  throughout  free  space. 

Ray  trajectories,  which  are  computed  by  the  simulation  and  shown  in 
7ig.  (24),  are  launched  from  the  emitter  at  0.02  deg  increments  between 
the  angles  of  *0.50  and  -0.5 0  deg  in  elevation.  This  ray  density  of  5° 
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rays/d eg  resulted  from  an  investigation  during  the  early  development  of 
the  simulation  which  indicated  that  the  field  strength  calculations,  while 
being  highly  sensitive  to  ray  spacings  greater  than  0.05  deg,  tended  to 
converge  rapidly  for  spacings  of  0.035  or  less.  A  ray  separation  of  0.02 
deg  was  then  selected  to  ensure  the  accuracy  and  consistency  of  the  simu¬ 
lation  field  strength  results. 

The  rays  in  Fig.  (24)  are  shown  with  an  artificial  upward  curvature 
that  results  from  plotting  the  earth’s  surface  along  a  linear  rather  than 
a  curved  axis.  This  curvature  is  actually  l/a  where  a  is  the  mean  radius 
of  the  earth.  If  the  rays  were  drawn  in  a  spherical  coordinate  system 
they  would  he  correctly  shown  to  be  straight  lines.  The  choice  of  a  rec¬ 
tangular  coordinate  system  was  made,  however,  for  ease  of  comparison  be¬ 
tween  simulation  results,  waveguide  calculations,  and  the  Guadalupe  Island 
measurements  of  Fig.  (20).  Another  point  to  be  made  is  that  the  rays  ap¬ 
pear  to  be  traveling  upward  at  a  rapid  rate  because  of  the  enormous  scale 
compression  along  the  abscissa  (approximately  22S  times  that  of  the  or¬ 
dinate  axis). 

?be  ray  shown  farthest  to  the  right  in  Fig.  (24)  defines  the  geomet¬ 
ric  horizon  seen  by  the  transmitter.  For  an  emitter  in  free  space  whose 
altitude  is  100  ft,  the  horizon  appears  0.18  deg  below  the  local  horizon¬ 
tal  at  the  transmitter,  given  that  the  earth's  surface  is  perfectly  smooth. 

-"•ays  launched  below  this  angle  will  be  reflected  from  the  earth  as  shown  in 
the  ray  trace  figure. 

Flgs-  (25)  and  (26)  show  the  refractlvity  profiles  and  ray  trajector- 
ies  for  the  same  emitter  using  the  standard  C.R.F.L.  exponential  atmos- 
.te-e  riven  by  Sqs.  (l?)  and  (129)*  A  comparison  of  Figs.  (24)  and  (26) 

**  °Ws  the  negative  N-gradient  of  the  exponential  atmosphere  (-13.65 


ATMOSPHERIC  REFRRCTIVI TY  PROFILE  REFRRCTIVI TY  ORflOIEHI  PROFILE 

FREE  SPACE  MODEL  FREE  SPRCE  MODEL 


Fig.  23.  1!  and  dN/dh  profiles  for  free  space 


TRUCE  FOR  THE  FREE  SPRCE  MODEL 


Fig,  24.  Ray  trace  for  a  100  ft  high  emitter  in  free  space 
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N-units/kft  at  a  height  of  100  ft)  causes  the  rays  to  refract  beyond  the 
geometric  horizon  such  that  the  distance  to  the  horizon  is  extended  ap¬ 
proximately  1.33  times  the  horizon  distance  in  free  space. 

Turning  now  to  the  Guadalupe  Island  duct,  Figs.  (2?)  through  (33) 
show  simulation  results  obtained  with  the  average  measured  refractivity 
profile  of  Table  (7)  and  Fig.  (21a).  Fig.  (2?)  shows  the  trajectories 
of  rays  launched  between  -+0.50  and  -0.50  deg  from  the  100  ft  high  emit¬ 
ter.  Only  those  rays  whose  elevation  angles  are  less  than  or  equal  to 
some  critical  angle  are  trapped  between  the  earth  and  the  1000  ft  top  of 
the  duct.  The  maximum  launch  angle  for  trapped  rays  may  be  found  from 
Snell's  law,  which  is  given  by  Eq.  (54)  for  spherically  concentric  media. 
Rewriting  Eq.  (54)  in  terms  of  ray  elevation  angles  (i.e.,  the  angular 
complement  of  3  and  fcQ)  gives 

rn(r)  cos  €  =  rQn(ro)cos  $o  (190) 


where  r,  n(r),  and  €  denote  values  at  the  point  of  refraction,  and  rQ, 
n(rQ),  and  €q  are  measured  at  some  known  point,  such  as  at  the  emitter. 
Setting  €  *  0  for  total  internal  reflection  (trapping)  at  the  maximum 
height  of  the  duct,  substituting  the  coordinate  change  of  Eq.  (60a)  and 
solving  for  €  ,  the  critical  launch  angle  for  trapped  rays  becomes 


where  h  and  hQ  are  the  altitudes  of  the  top  of  the  duct  and  the  transmit¬ 
ter  (measured  from  mean  sea  level)  and  a  ■  2.0925*10^  ft  is  the  mean  ra¬ 
dius  of  the  earth.  For  a  100  ft  transmitter  in  the  Guadalupe  Island  duct, 
the  maximum  launch  angles  are  |?o|  *  0.3&9  deg  £°r  ray  trapping  at  the 
top  of  the  duct  (h  =  1000  ft)  and  |€o|  =  0.255  deg  for  ray  trapping  at 
the  intermediate  800  ft  layer  in  the  duct.  Fig.  (27)  shows  that  ray  trap¬ 
ping  occurs  mainly  between  these  two  heights.  On  the  other  hand,  Eq. 


RAY  TRACE  FOR  THE  OUAOALUPE  ISLAND  DUCT 


Hay  trace  for  a  100  ft  emitter  in  the  Guadalupe  Island  duct 
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(191 )  indicates  that  ray  trapping  does  not  occur  at  the  2000  ft  high  layer 
shown  in  Figs.  (21a)  and  (22a)  for  any  value  of  €  . 

Fig.  (27)  shows  a  number  of  regions  where  there  is  a  distinct  absence 
of  rays,  particularly  the  area  above  the  duct  and  the  low  altitude  hole 
extending  from  10  to  50  naut  mi  within  the  duct.  The  hole  represents  the 
region  beyond  the  earth's  horizon  where  rays  are  unable  to  penetrate  and 
is  frequently  referred  to  as  the  earth  shadow  region,  Khile  geometric 
optics  predicts  that  there  is  no  field  present  in  either  of  these  regions, 
the  experimental  results  and  waveguide  mode  theory  calculations  presented 
in  the  following  section  indicate  that  just  the  opposite  is  true.  This 
discrepancy  arises  from  the  fact  that  classical  geometric  optics,  because 
it  does  not  make  explicit  use  of  wavelength  and  phase,  is  not  able  to 
solve  for  diffracted  fields  or  fields  resulting  from  evanescent  and  leaky 
modes  which  are  often  present  in  atmospheric  ducts.1,  12 '  1^'  1^  A  fur¬ 
ther  discussion  of  these  mode  types  and  their  fields  will  be  given  in  the 
latter  part  of  this  chapter. 

The  wavelike  ray  structure  shown  near  the  top  of  the  duct  in  Fig. 

(27)  forms  several  caustics,  or  regions  where  the  ray  density  increases 
rapidly  and  then  falls  off  abruptly  to  zero  (Chapter  IV).  Caustics  must 
be  excluded  from  consideration  since  they  violate  the  requirement  of  Eq. 
(115)  that  the  relative  spacing  between  adjacent  rays,  and  hence  the  ray 
density  in  a  volume  of  space,  must  not  change  appreciably  over  a  wave¬ 
length  of  distance. 

Kays  that  intersect  the  earth  in  Fig.  (27 )  are  reflected  from  sur¬ 
faces  that  are  locally  plane  and  tangent  to  a  smooth  spherical  earth  at 
the  point  of  reflection  (Chapter  V).  Thus  ray  divergence  due  to  reflect¬ 
ion  from  a  curved  earth  is  provided  for  in  the  ray  trace  diagram.  Surface 
roughness  becomes  a  consideration  only  when  calculating  the  field  of  a 
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reflected  wave  by  means  of  Fqs.  (184)  and  (185). 

Figs.  (28)  and  (29)  show  an  expanded  view  of  the  Guadalupe  Island 
refractivity  and  refractivity  gradient  profiles  of  Figs.  (21a)  and  (21b) 
and  the  ray  trace  plot  of  Fig.  (27)  for  heights  of  0  to  1000  ft  within  the 
duct,  since  the  region  above  the  duct  is  not  modeled  correctly  by  geomet¬ 
ric  optics.  Figs.  (30)  through  (33)  present  height  gain  curves  computed 
by  the  simulation  at  65,  170,  520,  and  3300  MHz  for  the  100  ft  high  emit¬ 
ter,  using  the  Guadalupe  Island  refractivity  profile  and  a  moderate  sea 
surface  (Ah  -  1  .1  ft)  model.  Vertical  reference  axes  are  drawn  at  20 
naut  mi  intervals  to  represent  the  zero  db  gain  level  of  field  strength 
relative  to  free  space  values.  A  scale  for  measuring  relative  field 
strength  is  given  in  the  upper  right  hand  corner  of  each  plot. 

The  most  notable  features  of  the  height  gain  curves  are  the  peaks 
and  nulls  which  result  from  mutual  interference  between  the  refracted  and 
reflected  wavefronts  propagating  through  the  duct.  Unlike  standard  micro- 
wave  waveguides,  the  duct  has  a  height  which  is  several  orders  of  magni¬ 
tude  greater  than  a  wavelength  of  radiation,  and  thus  is  capable  of  sus¬ 
taining  a  vertical  standing  wave  pattern  with  numerous  minima  and  maxima. 
Furthermore,  since  the  medium  within  the  duct  is  not  homogeneous  and  since 
the  boundaries  of  the  duct  are  considerably  different  (the  earth  surface 
causes  an  amplitude  and  phase  change  in  reflected  fields  while  the  more 
amorphous  upper  "boundary"  of  the  duct  does  not)  the  interference  patterns 
are  not  symmetric  either  with  respect  to  height  or  to  the  zero  db  refer¬ 
ence  axes.  Finally,  the  number  of  peaks  and  nulls  in  the  height  gain 
curves  is  seen  to  increase  with  emitter  frequency.  This  is  in  agreement 
with  elementary  optics  theory  which  states  that  the  distance  between  the 
minima  and  maxima  of  interfering  fields  is  inversely  proportional  to  fre¬ 
quency.  Thus  more  peaks  and  nulls  appear  as  this  distance  decreases. 


Simulation  results  for  the  500  ft  transmitter  are  shown  in  Figs. 

(34)  through  (33).  In  this  case  the  rays  are  launched  between  -HD. 54  and 
-0.45  deg  to  permit  the  simulation  to  detect  the  caustics  at  the  top  ef 
the  duct  and  thus  exclude  them  from  the  field  strength  calculations  (at 
least  two  adjacent  rays  must  lie  above  a  caustic  to  enable  its  location 
by  the  computer  program). 

The  ray  trace  of  Fig.  (34)  is  considerably  different  from  that  of 
Fig.  (29)  with  many  of  the  rays  trapped  between  the  emitter  and  the  bot¬ 
tom  edge  of  the  ducting  layer.  From  Snell's  law  and  Fig.  (34)  it  may  be 
shown  that  rays  having  0  <  <  O.O987  deg  are  trapped  between  500  and 

630  ft,  while  those  for  which  |=0J  <  0.411  deg  are  trapped  at  800  ft 

and  those  launched  at  | €q  |  <  0.490  deg  are  trapped  at  1000  ft.  The 

difference  in  ray  trapping  at  the  two  emitter  heights  may  be  explained 
by  3q.  (191).  Assume  for  the  moment  that  the  ratio  n(h)/n(ho)  does  not 
change  over  some  range  of  altitude  h.  3q.  (191)  states  that  as  h  becomes 
larger,  €q  will  become  smaller  until  at  some  point  the  right  hand  side  of 
the  equation  exceeds  unity  and  no  rays  may  be  trapped.  Thus  the  greater 
the  separation  between  the  layer  and  emitter  altitudes  (i.e.,  h  and  hQ) 
the  smaller  the  launch  angle  must  become  to  remain  trapped.  Of  course 
changing  n(h)  significantly  will  alter  this  result,  but  in  the  case  of 
trapping  at  the  800  and  1000  ft  levels,  the  ratio  of  n(h)/n(lv  )  remains 
relatively  constant  as  compared  to  the  ratio  of  (a  +  h)/(a  +  hQ)  when 
solving  2q.  (191)  at  hQ  -  100  ft  and  hQ  -  500  ft.  Similarly  3q.  (191) 
indicates  that  no  ray  trapping  will  occur  along  the  bottom  edge  of  the 
duct  for  the  100  ft  emitter,  as  confirmed  by  Fig,  (29)* 

Figs.  (35)  through  (33)  show  simulated  height  gain  curves  for  the  500 
ft  emitter  above  the  same  moderately  rough  (sea  state  3)  surface.  As  with 
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the  100  ft  emitter,  the  peaks  and  nulls  In  the  height  gain  curves  become 


more  pronounced  at  higher  frequencies. 

Finally,  the  simulation  results  using  the  trlllnear  ground  based 
duct  model  of  Pappert  and  Good hart  are  presented  In  Figs.  (39)  through 

(49).  Hays  are  again  launched  from  +0.50  to  -0.50  deg  in  elevation  for 
the  100  ft  emitter  and  from  **0.54  to  -0.46  deg  for  the  500  ft  emitter. 

A  comparison  of  critical  angles  for  ray  trapping  within  the  Guadalupe 
Island  and  trlllnear  duct  models  is  given  in  Table  (9). 


TA3L2  9 


7JI*  lirtUTTU-i 
TTJSSED  PAYS 


Emitter 

Trapping 

Maximum  Launch 

Model 

Height 

Height 

Angle 

(ft) 

(ft)  „ 

(deg) 

Guadalupe  Island 

100 

1000 

+  0.369 

800 

+  0.255 

500 

1000 

+  0.490 

800 

+  0.411 

630 

+  0.099 

Trlllnear 

100 

1000 

+  0.381 

800 

+  0.082 

500 

1000 

+  0.500 

800 

+  0.335 

650 

+  0.093 

Fig.  28.  N  and  dN/ dh  profiles  for  the  Guadalupe  Island  duct  (expanded  from  Fig.  (21a)) 
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Fig.  30.  Height  gain  curves  at  65  MHz,  from  Figs.  (28)  and  (29 
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Figs.  (28)  and  (29) 


GUADALUPE  ISLAND  DUCT  ABOVE  A  MODERATE  SEA  (STATE  =3)  HO  =  100  FT 


Fig.  32.  Height  gain  curves  at  520  MHz,  from  Figs.  (28)  and  (29) 
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Fig.  33-  Height  gain  curves  at  3300  MHz,  from  Figs.  (28)  and  (29) 


ROY  TRACE  FOR  THE  GUADALUPE  ISLAND  OUCT 


Fig.  34.  Ray  trace  for  the  500  ft  emitter  in  the  Guadalupe  Island  duct 
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Fig.  35.  Height  gain  curves  at  65  MHz,  from  Figs.  (28)  and  (34) 
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Fig.  36.  Height  gain  curves  at  170  MHz,  from  Figs.  (28)  and  (34) 
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Fig.  37.  Height  gain  curves  at  520  MHz,  from  Figs.  (28)  and  (34) 


GUADALUPE  ISLAND  DUCT  ABOVE  A  MODERATE  SEA  (STATE  =  3)  HO  =  500  F,T 


Fig.  38.  Height  gain  curves  a*  3300  MHz,  from  Figs.  (28)  and  (34) 


II  and  dll/dh  nrofiles  for  the  trilinear  /ground  based  duct  (expanded  from  Fig.  (21b)) 
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GROUND  BASED  DUCT  ABOVE  A  MODERATE  SEA  (STATE  =3)  HO  =  100  FT 
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Fig.  42.  Height  gain  curves  at  170  MHz,  from  Figs.  (39)  and  (40) 


Fig.  43.  Height  gain  curves  at  520  MHz,  from  Figs.  (39)  and.  (40) 


Fig.  44.  Height  gain  curves  at  3300  MHz,  from  Figs.  (39)  and  (40 ) 


RRY  TRfiCE  FOR  THE  TRIUNERR  GROUNO  BftSEO  DUCT 
HO  =  500  FT  1-0.54  TO  -0.46  DEGREE 


Fig.  45.  Kay  trace  for  the  500  ft  emitter  in  the  trilinear  duct 


GROUND  BASED  DUCT  ABOVE  A  NODERATE  SEA  (STATE  =  3J  HO  =  500  FT 
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Fig.  46.  Height  gain  curves  at  65  MHz,  from  Figs.  (39)  and  (45) 


GROUND  BASED  DUCT  ABOVE  A  MODERATE  SEA  (STATE  =3)  HO  =  500  FT 
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Fig.  47.  Height  gain  curves  at  170  MHz,  from  Figs.  (39)  and  (45) 
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igs.  (39)  and  (45) 


GROUND  BASED  DUCT  ABOVE  A  MODERATE  SEA  (STATE  =3)  HO  =  500  FT 


Fie.  49.  Height  gain  curves  at  3300  MHz,  from  Figs.  (39)  and  (45) 


The  preceding  simulation  results  axe  compared  in  this  section  to  the 
experimental  and  calculated  height  gains  presented  by  Pappert  and  Good- 
hart  for  the  ground  based  Guadalupe  Island  duct.  The  experimental  data 
for  the  100  ft  high  emitter  are  from  Pig.  (20)  and  the  waveguide  mode 
theory  calculations  are  based  upon  the  trilinear  refractivity  profile  of 
Table  (8)  and  Pig.  (21b).  The  height  gain  curves  from  Pappert  and  Good- 
hart,  shown  in  Figs.  (50)  through  (59),  are  for  fixed  receiver  heights 
of  100  and  500  ft  at  a  distance  of  120  naut  mi  from  the  transmitter 
(measured  along  the  earth's  surface),  with  selected  results  given  at  a 
distance  of  60  naut  mi.  Applying  the  principle  of  reciprocity  for  a  lin¬ 
ear  isotropic  medium,  these  results  are  compared  to  the  simulation  height 
gain  curves  where  the  transmitter  and  receiver  have  been  interchanged. 

As  mentioned  previously,  all  calculations  assume  a  sea  surface  of  moder¬ 
ate  roughness  (Ah  equals  1  ft  and  1.1  ft  in  the  respective  waveguide  and  - 
geometric  optics  results)  over  the  entire  th  between  Jan  Diego  and 
Guadalupe  Island. 

Pig.  (50)  shows  measured  and  calculated  height  gain  profiles  at  65 
MHz  for  a  transmitter  to  receiver  separation  of  60  naut  mi  with  the  re¬ 
ceiver  (or  simulation  transmitter)  at  500  ft  above  mean  sea  level.  Shown 

are  the  experimental  and  calculated  results  obtained  from  waveguide  the- 

1  14 

ory,  plus  calculated  height  gains  for  the  earth  diffracted  field  ’  and 

14  41 

the  field  scattered  by  the  troposphere.  ’  The  diffracted  and  tropo- 
scatter  fields,  which  are  labeled  "normal”  in  the  figure  legends,  are 
nonducted  fields  which  are  present  when  the  transmitter  is  below  the  hori¬ 
zon  (in  this  case  at  a  height  less  than  2200  ft).  These  fields  will  be 
shown  to  be  consistently  weak,  usually  ranging  from  20  to  60  db  below  the 
ducted  field.  The  experimental,  waveguide  theory,  and  normal  height  gain 
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curves  are  all  repeated  in  Fins.  (50a)  and  (50b)  to  avoid  congestion  when 
clotting  the  simulated  height  gains  obtained  from  the  Guadalupe  island 
and  trilinear  refractivity  profiles,  shown  in  Figs.  (50a)  and  (50b;,  re¬ 
spectively.  The  simulation  results  shown  in  this  case  are  from  the 
height  gain  curves  at  oO  naut  ni  given  in  Figs.  (35)  and  (46). 

In  examining  the  results  of  Fig.  (50)  it  is  apparent  that  geometric 
optics  yields  a  height  gain  profile  that  is  much  more  irregular  than  the 
curves  obtained  e;rperimentally  or  from  waveguide  theory.  This  may  be 
due  to  the  fact  that  the  geometric  optics  simulation  calculates  only  the 
fields  that  are  specularly  reflected  at  the  sea  surface,  and  thus  omits 
diffusely  scattered  fields  which  may  tend  to  "fill  in"  the  height  gain 


profile.  A  second  explanation  may  be  that  the  specular  reflection  co¬ 
efficient  of  the  sea  is  in  error  because  of  incorrect  electrical  proper¬ 
ties  assigned  to  the  sea  surface  in  Table  {?,).  Another  possibility  may 
lie  with  the  conditions  of  Iqs.  (77)  and  (73)  given  in  Chapter  IV  which 
state  that  geometric  optics  becomes  a  better  propagation  model  as  fre¬ 
quency  increases.  A  final  reason  may  be  that  since  geometric  optics  does 
not  make  explicit  use  of  wavelength  and  phase,  it  is  incapable  of  pre¬ 
dicting  the  existence  of  certain  types  of  modes  which  are  accounted  for 
in  waveguide  mode  theory.  It  will  be  shown  that  the  latter  two  explana¬ 
tions  are  the  most  likely  since  the  geometric  optics  results  do  improve 
considerably  at  higher  frequencies  and  since  geometric  optics  is  unable 
to'  calculate  the  field  above  the  duct  which  is  due  to  the  presence  of 
"leaky"  nodes  in  this  region.  If  the  lack  of  a  diffusely  reflected  field 
were  a  suitable  reason,  then  the  simulation  error  v.-ould  be  expected  to 
increase  rather  than  decrease  at  higher  frequencies,  since  the  sea  would 
appear  to  be  a  rougher  surface  and  thus  a  more  diffuse  reflector  at 
shorter  wavelengths.  Also,  an  incorrect  choice  of  electrical  properties 
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for  the  sea  would  not  be  expected  to  produce  a  significant  error  in  the 
Fresnel  reflection  coefficient  of  £c.  (ISO)  for  horizontally  polarized 
fields.  Little  change  occurs  in  the  reflection  coefficient  over  the 
range  of  snail  grazing  angles  encountered  (C  to  0.359  deg)  even  for  ex¬ 
tremely  large  variations  in  surface  electrical  characteristics  as  evi¬ 
denced  by  Figs.  (12)  through  (l?)  in  Chapter  V. 

"hile  both  the  experimental  data  and  waveguide  calculations  indi¬ 
cate  that  a  relatively  constant  field  is  measured  at  the  receiver  when 
the  transmitter  is  above  the  duct,  the  reciprocal  height  gain  curve  ob¬ 
tained  from  geometric  optics  shows  no  field  present  for  this  condition. 
Waveguide  mode  theory  predicts  that  such  a  field  does  exist  which  is  the 

result  of  the  coupling  of  energy  from  the  transmitted  signal  into  the 

1  14 

duct  by  means  of  "leaky"  modes. *’  If  the  transmitter  and  receiver 
were  reversed,  as  they  are  in  the  simulation,  then  this  field  would  ex¬ 
ist  above  the  duct  because  of  energy  leakage  by  the  same  types  of  modes. 
However  as  the  ray  traces  of  Figs.  (29),  (39-),  (40),  and  (b-5)  show,  there 
is  nothing  in  the  geometric  optics  solution  to  suggest  energy  leakage 
along  the  top  of  the  duct, 

To  aid  in  describing  energy  leakage  through  a  ducting  layer,  consid¬ 
er  the  comparison  of  a  tropospheric  duct  to  a  dielectric  slab  waveguide. 

In  the  case  of  the  dielectric  waveguide,  the  field  may  be  resolved  into  a 
sun  of  elementary  waves  or  modes  which  axe  guided  along  the  slab  boundar¬ 
ies  with  little  or  no  attenuation  in  the  direction  of  propagation  and 
’■'ith  an  exponential  decay  in  the  direction  normal  to  the  outside  of  either 
soundary.  However  z  duct  has  no  well  defined  upper  boundary,  such  as  a 
discontinuity  in  the  refractive  index,  and  consequently  has  propagation 
characteristics  which  are  different  from  those  of  a  simple  dielectric 
waveguide .  Hode  theory  predicts  that  two  types  of  nodes,  commonly  referred 
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to  as  "trapped''  and  "leaky"  modes,  may  exist  within  a  meteorological 
duct.  Trapped  modes  are  unattenuated  modes  of  propagation  which  are 


analogous  to  those  of  the  dielectric  slab.  Leaky  modes,  while  guided 
within  the  duct  along  the  earth's  surface,  are  allowed  to  propagate  ob¬ 
liquely  to  the  duct  in  the  region  above  the  superrefractive  layer.  En¬ 


ergy  is  then  coupled  into  or  out  of  a  duct  when  one  or  more  leaky  modes 

T  i  o  T  - 

is  strongly  excited  by  a  nearby  transmitter.- ’ 


figs.  (51)  and  (52)  show  results  at  65  II is  for  a  120  naut  mi  trans¬ 


mitter  to  receiver  separation  with  the  receiver  at  100  and  500  ft 
heights.  Again  Figs.  (51a)  and  (52a)  give  simulated  height  gain  curves 
using  the  Guadalupe  Island  refractivity  profile,  and  Figs.  (51b)  ana 
(52b)  show  simulation  results  obtained  with  the  trilinear  duct  profile. 
As  before,  the  geometric  optics  height  gain  curves  are  highly  irregular, 
especially  those  using  the  trilinear  refractivity  profile.  Also,  the 
simulated  height  gains  shown  in  these  figures  are  generally  5  to  20  db 
greater  than  those  of  either  the  experimental  or  waveguide  curves.  It 
is  quite  evident  that  waveguide  mode  theory  has  a  far  better  agreement 
with  the  experimental  measurements  at  65  .Hz  than  does  geometric  optics. 

Fappert  and  C-oodhart  indicate  that  equipment  calibration  errors 
have  sometimes  affected  the  measured  data,  resulting  in  this  instance 
in  a  larger  than  expected  difference  between  the  experimental  and  mode 
theory  height  gain  curves  shown  in  Fig.  (52).  Fvidence  of  equipment 
error  appears  in  Figs.  (51)  and  (52),  where  the  measured  field  at  500 
ft  in  Fig.  (51)  is  not  the  same  as  the  field  at  100  ft  in  Fig.  (52), 
thus  indicating  a  violation  of  reciprocity. 

2igs.  (53)  and  (5*0  show  experimental  and  calculated  height  gain 
curves  at  1?0  Ills  for  a  distance  of  12C  naut  mi  with  the  receiver  height 
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at  100  and  500  ft.  The  simulation  results  are  still  in  much  poorer  agree¬ 
ment  with  the  measured  data  than  are  the  waveguide  calculations,  however 
the  geometric  optics  results  of  Fig.  (53)  do  show  a  two  mode  interference 
structure  which  is  slightly  displaced  from  th°  same  two  mode  pattern 
given  by  waveguide  theory.  Pappert  and  Goodhart  note  that  the  number  of 
modes  required  in  the  waveguide  calculations  ranges  from  a  single  mode 
at  65  MHz  and  two  modes  at  170  MHz,  to  nearly  a  hundred  modes  at  3300 
MHz.  They  also  explain  that  the  sharp  null  at  1200  ft  in  the  experimen¬ 
tal  data  of  Fig.  (53)  is  most  likely  the  result  of  a  transmitter  and  re¬ 
ceiver  antenna  misalignment,  and  is  therefore  not  to  be  interpreted  as  a 
real  null. 

Figs.  (55)  and  (56)  show  height  gain  curves  at  520  MHz  for  a  range 
of  120  naut  mi  with  receiver  heights  of  100  and  500  ft.  The  simulation 
results  give  field  strengths  that  are  within  5  to  10  db  of  the  experimen¬ 
tal  and  mode  theory  results,  which  is  a  considerable  improvement  over  the 
previous  cases  at  the  lower  frequencies.  Agreement  is  generally  best  be¬ 
tween  the  simulation  results  using  the  Guadalupe  Island  profile  and  mode 
theory  results,  although  as  Fig.  (55a)  shows,  some  of  the  interference 
lobes  do  not  appear  in  the  geometric  optics  height  gain  curve  between  the 
altitudes  of  400  and  600  ft.  The  deep  null  at  1100  ft  is  again  consid¬ 
ered  to  be  the  result  of  antenna  misalignment. 

Fig.  (56)  shows  waveguide  calculations  for  cases  where  the  modal 
equation  is  solved  at  the  ground  (l  =  0)  and  at  the  bottom  edge  of  the 
ducting  layer  which  is  at  a  height  of  600  ft  (D  =  600).  Pappert  and 
Goodhart  state  that  at  this  frequency  modes  exist  which  are  either  earth 
detached  or  evanescent  at  the  ground.  The  height  gain  curve  for  D  =  600 
includes  these  additional  modes  and  is  in  better  agreement  with  experi¬ 
mental  results  by  eliminating  the  fine  structure  of  the  D  ■  0  curve. 
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"’hile  the  .-node  calculations  in  Pic.  (55)  are  for  D  =  600,  deep  nulls  ap¬ 
pear  in  the  waveguide  height  gain  curve  since  the  additional  inodes  are 
sufficiently  evanescent  at  the  100  ft  receiver  height  so  as  not  to  af¬ 
fect  the  node  sum  by  any  appreciable  amount. 

Since  geometric  optics  solves  for  the  direction  of  energy  propaga¬ 
tion  through  a  medium,  it  is  unable  to  model  the  distribution  of  energy 
contained  in  the  nonpropagating  evanescent  modes  existing  within  a  duct. 
However  ray  optics  is  capable  of  describing  earth  detached  propagation 
as  evidenced  by  the  ray  traces  of  Figs.  (29),  (34),  (40),  and  (45). 

These  figures  all  show  a  number  of  rays  which  travel  below  the  ducting 
layer  without  Intersecting  the  earth's  surface.  Such  elevated  rays  rep¬ 
resent  earth  detached  propagation  which  begins  to  appear  in  the  simulated 
height  gain  curves  at  higher  frequencies.  In  the  520  MHz  case,  geometric 
optics  produces  results  which  lie  between  the  mode  theory  height  gains 
for  D  *  0  and  D  -  600.  It  should  be  noted  that  neither  ray  optics  nor 
mode  theory  gives  results  that  are  in  very  good  agreement  with  the  experi¬ 
mental  data  in  Pigs.  (55)  and  (56) .  This  discrepancy,  which  is  mentioned 
by  Pappert  and  Goodhart  for  the  waveguide  calculations,  is  suspected  to  be 
due  to  temporal  or  spatial  fluctuations  in  the  duct  layer  at  the  time  of 
measurement. 

Fig.  (57)  shows  results  at  330 0  MHz  for  a  range  of  60  naut  mi  and  a 
receiver  height  of  100  ft.  The  pair  of  solid  curves  for  the  experimental 
data  represent  the  envelope  of  measured  field  strength  at  this  range.  The 
results  obtained  from  geometric  optics  are  in  very  good  agreement  with 
both  the  experimental  and  waveguide  theory  height  gain  curves,  mode  the¬ 
ory  predicts  a  large  decrease  in  the  field  above  300  ft  because  of  the 
destructive  interference  in  this  region  of  nearly  100  modes.  Pappert  and 
Goodhart  point  out  that  this  phasing  would  be  eliminated,  however,  by  a 
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turbulence  in  the  ducting  layer  of  the  trilinear  refractivity  model. 

Pigs.  (58)  and  (59)  show  results  at  3300  MHz  for  a  transmitter  to 
receiver  separation  of  120  naut  mi  with  the  receiver  at  100  and  500  ft. 
'..'hile  the  simulation  results  fall  between  the  experimental  and  waveguide 
height  gain  curves  for  the  100  ft  receiver  of  Pig.  (53),  both  geometric 
optics  and  mode  theory  predict  fields  that  are  5  to  10  db  higher  than 
the  field  measured  at  this  altitude.  Unless  both  mode  theory  and  geomet¬ 
ric  optics  are  in  error,  the  most  likely  explanations  are  that  either  an 
antenna  alignment  error  still  continued  to  exist  or  the  duct  did  not  re¬ 
main  stationary  and  horizontally  homogeneous  at  the  time  of  measurement. 

In  the  case  of  the  500  ft  high  receiver  the  simulated  height  gain  curves 
given  for  the  Guadalupe  Island  profile  in  Fig.  (59a)  again  appear  to  be 
in  better  agreement  with  the  mode  theory  calculations  than  those  simu¬ 
lated  with  the  trilinear  profile  in  Fig.  (59b).  The  experimental  results, 
while  lacking  the  fine  structure  of  either  the  ray  optics  or  waveguide 
height  gain  curves,  are  shown  to  be  an  approximate  average  of  the  calcu¬ 
lated  curves  in  Fig.  (59).  This  may  be  due  to  either  a  change  in  the 
duct  refractivity  which  is  not  accounted  for  in  either  duct  profile  model 
or  the  existence  of  other  earth  detached  or  evanescent  modes  which  were 
not  included  in  either  the  mode  theory  or  geometric  optics  results. 
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Fig.  52.  Height  gain  curves  at  65  MHz,  for  h  =  500  ft  and  x  =  120  naut  mi 
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.  54,  Height  gain  curves  at  170  i'iHz,  for  h  -  500  ft  and  x  =  120  naut  mi 
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Fig.  55*  Height  rain  curves  at  520  MHz,  for  h  =  100  ft  and  x  =  120  naut  mi 
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CHAPTH.  VII 

The  results  presented  in  this  paper  indicate  that  classical  geomet¬ 
ric  ortics  yields  moderately  fair  predictions  of  the  strength  of  a  trop- 
ospherically  ducted  field  at  frequencies  below  100  to  200  KHz,  and  rea¬ 
sonably  good  estimates  at  frequencies  above  500  KHz.  The  fields  calcu¬ 
lated  from  geometric  optics  were  generally  higher  and  more  irregular 
than  those  obtained  experimentally  or  from  waveguide  mode  theory  at  65 
and  170  iris,  although  at  170  KHz  both  mode  and  ray  theories  differed  at 
times  from  the  measured  field  by  as  much  as  20  db.  These  errors  may  be 
related  to  problems  in  the  alignment  and  calibration  of  the  measurement 
equipment  or  in  temporal  and  spatial  fluctuations  in  the  ducting  layer. 

In  the  case  of  the  geometric  optics  results  a  further  explanation  may 
lie  with  the  restrictions  of  Tas.  (77)  and  (7S)  in  Chapter  IV,  which 
state  that  ray  theory  is  useful  as  a  propagation  model  only  at  high  fre¬ 
quencies.  Nonetheless,  ray  optics  appears  to  provide  an  order  of  magni¬ 
tude  approximation  to  the  field  as  low  as  65  KHz ,  which  could  conceivably 
be  extended  down  to  3 0  KHz  where  atmospheric  refraction  becomes  the  pri¬ 
mary  mechanism  for  long  distance  radio  propagation. 

3y  contrast,  the  fields  calculated  from  geometric  optics  at  320  and 
3300  KHz  were  generally  within  5  to  10  db  of  both  the  experimental  and 
waveguide  theory  results.  In  some  of  the  examples  shown  at  these  higher 
frequencies,  ray  theory  was  in  better  agreement  with  waveguide  theory  than 
’•’ith  the  experimental  data.  Again  this  difference  may  have  resulted  from 
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equipment  error  or  inhomogeneities  in  the  duct  layer.  Finally,  the  ray 
theory  calculations  which  had  the  closest  agreement  with  either  experi¬ 
mental  or  mode  theory  results  were  those  derived  from  the  Guadalupe  Is¬ 
land  refractivity  profile.  The  only  plausible  reason  for  this  outcome 
is  that  the  Guadalupe  Island  profile,  by  including  several  variations  in 
the  refractivity  structure  both  within  the  ducting  layer  and  near  the 
earth  surface  which  were  not  modeled  by  the  trilinear  profile,  simply 
was  more  representative  of  the  Guadalupe  Island  duct. 

Limitations  that  are  inherent  to  classical  ray  theory  have  been 
noted  throughout  this  paper.  These  include  the  inability  to  model  dif¬ 
fracted  fields  and  certain  propagating  (leaky)  and  nonpropagating  (eva¬ 
nescent)  modes  which  are  present  in  a  tropospheric  duct.  Regions  in 
which  these  effects  are  of  concern,  such  as  the  earth  shadow  region  and 
the  area  above  a  duct,  must  then  be  treated  by  the  diffraction  and  wave¬ 
guide  mode  theories  of  physical  optics. 

Despite  its  obvious  limitations,  geometric  optics  provides  a  useful 
qualitative  as  well  as  quantitative  description  of  the  effects  of  atmos¬ 
pheric  refraction  on  wave  propagation.  Given  a  distribution  of  atmos¬ 
pheric  refractive  index,  the  ray  trajectories  of  geometric  optics  yield 
an  easily  understood  representation  of  wavefront  propagation  through  lay¬ 
ered  atmospheres  with  a  minimum  of  computational  effort.  Furthermore,  un¬ 
like  waveguide  mode  theory  where  large  numbers  of  modes  may  be  required  to 
obtain  a  field  solution,  geometric  optics  is  able  to  make  rapid  and  effi¬ 
cient  field  calculations  for  different  frequencies,  polarizations,  and 
surface  properties  based  upon  a  single  set  of  ray  trajectories.  This  a- 
bility  to  lend  a  physical  interpretation  to  refraction  effects  plus  the 
®ase  and  efficiency  of  computation,  make  ray  theory  an  attractive  alter¬ 
native  to  the  more  laborious  mode  theory  of  propagation,  especially  when 


several  layers  are  present  in  the  atmosphere.  Moreover,  the  results  pre¬ 
sented  in  this  thesis  indicate  that  geometric  optics,  when  properly  un¬ 
derstood  and  applied,  offers  a  reasonably  accurate  method  for  determining 
the  strength  of  a  field  in  an  inhomogeneous  atmospheric  structure,  such 
as  a  tropospheric  duct,  which  is  comparable  in  realism  to  that  of  wave¬ 
guide  mode  theory. 


This  Appendix  provides  a  description  plus  instructions  for  use  of 
the  atmospheric  refractivity  computer  simulation  ( ATREF )  discussed  in 
Chapter  V.  Program  ATREF  is  written  in  Extended  FCRTHAL'  language  for 
the  CJC  6600  digital  computer  system  and  requires  approximately  65,000 
octal  words  of  memory.  Approximately  12  to  15  seconds  of  central  proc¬ 
essor  time  are  required  to  calculate  the  ray  trajectories  and  height 
gains  for  an  emitter  at  one  frequency.  Program  input  is  given  in  stan¬ 
dard  card  image  form  and  output  is  provided  in  both  printed  and  on-line 
CALC 0? 2?  plotted  formats,  depending  on  user  selection.  A  description  of 
the  executive  program  and  seven  subroutines  of  the  computer  simulation 
is  given  below. 

ATREF  is  the  main  executive  program  which  controls  the  input/output 
functions,  refractivity  model  selection,  and  calculation  of  the  ray  tra¬ 
jectories  and  height  gain  curves.  Refractivity  profiles  may  consist  of 
either  the  stored  free  space  or  exponential  CIIPL  models  given  by  Sas. 
(l6),  (129),  and  (130),  or  an  input  profile.  Input  profiles  are  stored 
as  piecewise  linear  functions  connecting  the  data  points  of  the  profiles. 

R2FRCT  performs  the  printing  and  plotting  of  the  refractivity  and 
refractivity  gradient  profiles, 

RrXrrc  sets  up  the  integration  of  the  ray  trajectory  and  time  of 
propagation  differential  equations  given  by  3qs.  (131)  through  (133). 

'•'hen  a  ray  crosses  a  boundary  between  two  pieceuise  linear  segments  of 
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an  input  refractivity  profile  or  when  a  ray  intersects  the  earth,  REIKTG 
integrates  to  the  boundary  by  means  of  a  variable  step  size  interpola¬ 
tion  algorithm.  The  ray  equations  are  then  reinitialised  at  the  boundary 
and  integration  continues  until  another  boundary  is  crossed,  repeating 
the  interpolation  process. 

RK  integrates  the  ray  equations  by  means  of  a  fourth  order  Runge- 
Kutta  algorithm  using  either  an  integration  step  size  of  D3L»XD3LTA/lC, 
where  Fd/RLTn  is  the  distance  interval  along  the  earth's  surface  for 
printing  and  plotting  the  height  gain  curves,  or  a  step  size  computed 
by  the  interpolation  operation  in  subroutine  RKIKTG . 

ATI-133  computes  atmospheric  refractivit'*  at  the  ray  height  and  sets 
up  the  ray  equations  for  use  in  subroutine  RK, 

3CATT  calculates  the  complex  Fresnel  reflection  coefficient  and 
surface  roughness  factor  for  use  in  Rqs.  (153) ,  (159),  and  (1B7). 

HCL'iXN  computes  the  relative  field  strength  and  power  density  by 
means  of  K'as.  (137)  through  (loO)  in  Chapter  V.  The  height  gain  is  ob¬ 
tained  from  field  calculations  at  100  overlapping  "window"  or  altitude 
increments  (of  height  dly  in  Chapter  V)  between  the  highest  and  lowest 
rays,  legions  that  are  not  correctly  modeled  by  geometric  optics,  such 
as  earth  shadow  regions  and  caustics,  are  omitted  in  the  height  gain 
plots  and  represented  in  the  printed  output  by  asterisks  (*). 

PLOTTR  draws  and  labels  all  plot  axes  and  prints  two  lines  for  lab¬ 
eling  above  each  plot.  PLOTTR  also  prints  a  one  line  banner  preceding 
the  plots  for  additional  descriptive  purposes  (e.g. ,  the  date),  and  pro¬ 
vides  for  scaling  the  plots  from  their  normal  sizes  (4x5  inches  for  X 
and  dIC/dh  profiles  plots  and  10x5  inches  for  ray  traces  and  height  gain 
curves). 


on  oo  onooooonnoonoooooa 


141 


PWJWH  ATREF  (INPUT ,  OUTPUT  ,  PLOT ) 

- —THIS  PROGRAM  IS  A  GEOMETRIC  OPTICS  MODEL  OF  HAVE  PROPAGATION 

- THROUGH  AN  INHOMOGENEOUS  ATMOSPHERE  HAVING  A  VERTICALLY 

- STRATIFIED  INO^X  OF  REFRACTION.  THE  “ROGRAH  CALCULATES  THE 

— — OIRECTION  OF  WAVEFRONT  PROPAGATION  BY  SOLVING  THE  EJLER- 
— — LAGRANGE  EQUATIONS  OF  RAYS  NORMAL  TO  INCREMENT6.  WAVEFRONT 

- SURFACES.  THE  RAY  TRAJECTORIES  ARE  THEN  USED  TO  COMPUTE 

- THE  RELATIVE  EMITTER  FIELD  STRENGTH  OR  POWER  OcMSITY 

- (NORMALIZED  TO  FREE  SPACE)  AS  A  FUNCTION  OF  ALTITUOE  AND 

- OISTANCE  ALONG  THE  EARTH'S  SURFACE.  FIELDS  WHICH  ARE 

'——REFLECTED  FROM  THE  EARTH  ARE  ATTENUATED  BY  A  FRESNE. 

- REFLECTION  COEFFICIENT  AND  A  SURFACE  ROUGHNESS  r ACTOR. 

— - THE  ELEVATION  ANGLE  AND  TIME  OF  PROPAGATION  ARE  CALCULATED 

- ALONG  EACH  RAY  PATH  TO  DETERMINE  THE  DIRECTION  OF  THE 

— — WAVEFRONT  PROPAGATION  VECTOR  AN3  THE  PHASE  RELATIONSHIP 
——BETWEEN  INTERFERING  WAVEFRONTS  rOR  THE  FIELD  STRENGTH  AND 
- POWER  DENSITY  COMPUTATIONS. 

COMMON  /ONECOM/  CRH, CRX, CRG, OTR, JAREA , NO ATA ,C, Ct ,C2 

1  ,RHO(51.30)  ,PHI(51,3!))  , PI , 3MF, CNAUT . AHS 

2  ,*RQ,TPW,NHV,NSL,NRMS,A8SRH, PHASE, THETA 

COMMON  /TKOCOM/  RA, 3ELX, AHO.ELOt 5i> , JRA Y , EXNAX , AHMAX.XFINAL 
COMMON  /THRCOH/  NELO, JPLOT ,IPLAC£«VX(3G) . JCASE 

1  • H (SI  ,30) ,Gt51«30)*HN(5G)«.RN(50),S(51,3G) 

2  , KSIZ.YSIZ, JPLT.XMIN, YMIN, XDI V, YDIV 

3  , XNAX, YMAX ,MR£F , NR AY , NPRO , NPLOT , NCR AD 

A  ,<REF,tCGRAO,KP.OT 

DIMENSION  P0L(2I,TER(4), REFMOO (2 , 3 ) 

DATA  P3L/10HHORIZONT AL,10HVERTICAL  / 

OATA  TER/10HSEA  WATER  ,10HORV  GROUND, 1  DMA VG  GROJND, 10HWET  GROUND/ 
OAT A  REFMOO/10HFREE  SPACE, 10H  MODEL  , 13HEXPONENTI A, 

1  10HL  MODEL  , 10HIN*UT  PROF.lOHlLE  MODE.  / 

OATA  DTR.RA,CM,PI/1.7«i53E-02,2.3925E*07,2.998EO8,3.14159/ 

OATA  CMF , CNAULT /3. 281, 6076.0/ 

SET  UP  THE  NUMBER  0-  PROGRAM  RUNS 
READ  30C.NCASE 
00  200  ICASE*1«NCAS£ 

READ  AMO  PRINT  THE  INPUT  OATA 
READ  300, NPRO 
READ  30C ,NOATA 

READ  310,1  CHN(IOI,  RN(IOM,  10*1,  MDATA) 

READ  320, AHS, AHO, AH MAX 

READ  3?0,XDELTA«XFINAL 

READ  320 ,EL0S1 ,ELOS2 

READ  320 ,FRQ, PW 

READ  300,NHV,NSL,MRMS 

READ  330,RREF,KGRAD,KRAY,KPLOT 

READ  3)9, NREF,NGRAO,NRAY, NPLOT 

PRINT  330 

.  PRINT  340, NCASE, NPRO, (REFMOO ( IRM, NPRO*l ) ,IRM=1,2I 
PRIM”  350 , AHS, AHO, AHMAX 
PRINT  3SO»XOELTA,XFINAL 
PRINT  3TO,ELOS1,ELOS2 

PRINT  ?SO,FRQ,PW,POw(NHV>,TE.RtNSL),NR.MS 


Fig.  61.  Program  listing 
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PRINT  J90  »  KR£F,KGRAO, KRAY , KPLOT 
PRINT  450,NREF,NGRA;),NRAY,NPLOT 
C 

C  CALCC-ATE  ZERO  MEAN  SEA  LEVEL  REFRACTI VITY  IF  NOT  ALREAOY  KNOWN 

NTEHP=NDATA 
IF  (HN( NOATA) )  6,6,5 

6  SLOP£=TRN(NOATA-l»-RNtNOATAM  /  HN INDAT  A-l) -HN  (NO  ATA)  ) 

RNINOAT  A-*-l )  =RN< NOATA)  -SLOP£*HN (NOATA) 

HNINOATA-H)  =0.0 
NTENP=NDATA*1 
C 

C  SET  UP  INITIAL  ATMOSPHERIC  REFRACTIVITY  CONSTANTS 
6  IF  (NPRO-l)  10,15,23 
C 

C  FREE  S=»AC£  MOOEL 
10  51*0.0 
C2=D. 0 
60  TO  25 
C 

C  EXPONENTIAL  MOOEL 
15  Cl=313.0 

C2=0. OOGG4366  - 
GO  TO  25 
C 

C  PIECE-WISE  LINEAR  MOOEL 
20  C1=RN<NTEHPI 

C2MAL0GIRN  INTEMP) ) -ALOS IRNtl ) )  1  /HNI1) 

C 

C  SET  UP  INITIAL  CONDITIONS 
25  C=CHF*CM 
NELO=50 

CNELO=( ELOS1-ELOS2) *FLOAT (NELO) 

OENOT*l. /FLOAT (NELO) 

DELX=CNAUT*XDELTA/13. 

EXMAX=CHAUT*XFINAL 

FR3=10)0000.0*FRQ 

PL4CE=EXHAX/OELX 

TPW=3H/100000  0.0 

ICALC*<PLOT*NRAY*NP_OT 

IPLACE*IFIX|PLACE)»l 

IPL0T*NREF«-NGRA0*NRAY+NPL0T 

JCAS£*ICASE 

JPLT=0 

NELO=IcIXtCNELO)-*l 

c 

C  CALL  ROUTINES  IF  THERE  IS  A  PRINTOUT  OR  PLOT  OF  THE 
C  REFR4GTI VITY  PROFILE 
JPLOT*l 

IF  (NR2F)  30,30,31 

30  IF  IKREF).  3 3,33,32 

31  CALL  PLOTTR 

32  CALL  REFRCT 
C 

C  CALL  ROUTINES  IF  THERE  IS  A  PRINTOUT  OR  PLOT  OF  THE 
C  REFRACTIVITY  GRADIENT 

33  J»LOT«2 

IF  (NGTAO)  34,34,35 


Fig.  6l.  (continued) 
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34  IF  (KGRAOJ  37,37,36 

35  CALL  P.OTTR 

36  CALL  REFRCT 

MiK  I-  THERE  ARE  ANY  FURTHER  CALCULATIONS 

37  IF  (IC4LCI  200,260,00 

CHECK  1-  THERE  IS  T3  BE  A  PLOT  OF  THE  RAY  TRACES.  IF  SO, 

CALL  ROUTINE  TO  SET  UP  THE  PLOT  AXES 
40  IF  (NRAYJ  44,00,42 
42  JPL0T=3 

CALL  PLOTTR 
44  CONTINJE 

SET  U°  A  LOOP  TO  CALCULATE  THE  ALTITUDE  PROFILE  OF  EACH  RAY 
DO  100  1  =  1, NELO 

INITIALIZE  FOR  INTE3RATIOM 

CRH«AH0 

CRKCQ.3 

ICCC-I *1 

ELOU)  =ELOSl-DENOT*fLOAT(ICCC» 

EANS=DTR»ELO(U 

CRG=( (RA*CRH)/RAI*SIN(EANS»/COS(EANG> 

JRAY=I 

CALL  ROUTINE  TO  COMMUTE  RAY  TRACES  ANO  PROPAGATION  TIMES 
CALL  RKINTG 
100  CONTINUE 

RESET  PARAMETER  VALJES 

OELX*13.*OELX 

IPLACE=I°LACE/10*1 

CALL  ROUTINES- IF  THERE  IS  A  PRINTOUT  OR  PLOT  OF  THE 
RELATIVE  FIELD  STRENGTH  OR  RELATIVE  POWER  DENSITY 
120  JPLOTs* 

IF  (NPLOT)  130,130,140 
130  IF  (KP.OTI  200,200,150 
140  CALL  PLOTTR 
150  CALL  H3AIN 
200  CONTINUE 

CALL  LIBRARY  ROUTINE  FOR  ON-LINE  PLOTTING 
IF  (IPLOT)  220,220,210 
210  CALL  P.OTE 
PRINT  430 
220  CONTINJE 

300  FORMAT  (6151 
310  FORMAT  (2F15.7I 
320  FORMAT  (5F15.7I 

330  FORMAT  ( 1H1 ,2X ,* ATMOSPHERIC  RADIO  REFRACTIVITY  COMPUTATIONS*// > 
340  FORMAT  (2X,*  NCASE=*,I6/ 

1  2X, *  NPR0**,I6/ 

2  2X , *  NODE***, 5X, 2A10) 

350  FORMAT  (2X,*  AH3**,F10.2,4X,*FEET*/ 
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1  2X,*  AH3=»,F10.2,4X,*FE£T*/ 

2  2X  »*  AHMAX=*,F10.2,4X,*FE£T*> 

360  FORMAT  (ZX  t*  XDELTA=*,F10.2,4X,*NAUT  MI*/ 

1  2X»  *  XFINAU**,F10.2,4X,*NAUT  HI*> 

370  FORMAT  t2X,*  ELOSl=*, F10 . 2, 4X, *  DEG*/ 

1  2X,*  ELOS2=*,F10.2,4X,*OEG*» 

300  FORMAT  C2X,*  FR3=*» F10 . 2, 4X, *MHZ*/ 

1  2X  »  *  PW=*»F10.2,4X.*MICROSEC*/ 

2  2X i *  POLAR**, 5X, AID/ 

3  2X *  *  TERRAIM**,5X,A1J/ 

4  2X,*  NRM3=*,I6> 

390  FORMAT  (2X,*  KREr=*,I6/ 

1  2X«  *  KGRAD=*« 16/ 

2  2X,*  KRA  f  =  * , IS / 

3  2X » *  KPLOT  =** 16 ) 

400  FORMAT  (2X,*  NREF=*,I6/ 

1  2X,*  N6RA3=*,IS/ 

2  2X»*  NRAT=*,I6/ 

3  2X,  *  NPLOT=*,I6> 

430  FORMAT  <5X,*ENO  OF  FILE  uN  PLOTTER  TAPE*! 
END 
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SU3R0UTINE  REFRCT 
C 

C~ - THIS  ROUTINE  C ALCUL A  TES  AND  PLOTS  THE  REFRACTIVITY  AND 

C - REFRACTIVITY  GRAOIENT  PROFILES 

C 

CONNON  /ONECOM/  CRH,  CRX , CRG.OTR,  JA REA,  NO  A7.\,C, Cl  ,02 

1  «  RHO (51,30 , PHI <51,30 ) , PI , CMF , CN AUT , AHS 

2  ,cRD,TPw,NHV,N3L,NRHS,ASSRH, PHASE, THETA 

COHN ON  /TWOCOM/  RA , 3ELX , A HO,£LO < 5 1 ) , JR  AT , EXHAX , A HM AX , XFINAL 
COMNON  /THRCOM/  NELO, JPLOT , IPLACE , VX < 30 > . JCASE 

1  »H(51,30),G(51*3GJ»HN<53),RN(5D)»S(51*30) 

2  ,XSIZ*YSIZ,JPLT,XMIN»YHIN,XDIV»YDIV 

3  ,<NAX,YMAX,NREr, NRA Y , NPRO, NPLOT , NGRAD 

4  ,  <REF ,  KGRAO,  KPi-OT 
DIMENSION  X  (150) , Y (150  I 

C 

C  PRINT  HEADING  IF  THERE  IS  A  PRINTOUT 
IF  CJPlOT-1)  10,10,14 
10  IF  (KREFI  18,18.,12 
12  PRINT  110 
GO  TO  IS 

14  IF  (KGRAO)  13,18,16 
16  PRINT  120 

C 

C  SET  UP  INITIAL  CONDITIONS  TO  CALCULATE  REFRACTIVITY  AND 

C  .  REFRACTIVITY  GRADIENT  VERSUS  ALTITUDE 

15  IF  (NPRO-1)  20,20,33 
20  10=51 

DILH=AHNAX/50. 

A=AHMAX*DELH 
GO  TO  40 
30  IO=NOATA 
L  =  0 
M*0 

IF  (HN (1 ) - AHMAX)  35,60,60 
35  OELH=(6HHAX-HN(ll)/25. 

A=AHNAX*OELH 

10=24 

N=I0-1 

C 

C  SET  JP  A  LOOP  TO  PRINT  AND  PLOT  THE  FREE  SPACE  AND  EXPONENTIAL  PROFILES 
40  CONTINUE 

DO  55  I=1,ID 
A=A-DELH 

XCI)*Cl*EXP(-Ci*A) 

YCI)=A 

IF  (JPLOT-1)  42,42,46 
C 

C  PRINT  THE  REFRACTIVITY  PROFILE 

42  X1=X(I)  - 

T1=YCI) 

IF  (KREF)  44,44,43 

43  PRINT  140,1, Yl, XI 

44  IF  (N-REF)  54,54,50 
C 

C  PRINT  THE  REFRACTIVITY  GRAOIENT  PROFILE 

46  IF  CI-IOJ  47,54,54 
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47  Xl=1003.*(X(I)-X(I+l)l/(Y<I)-Y(I+l)) 

Yl=IYm+V(I+lM/2 .3 

IF  IKGRAD)  49,49,48 

48  PRINT  140, 1, Yl, XI 

49  IF  f NGR AD)  54,54,50 

PLOT  THE  REFRACTIVITY  AND  REFRA2TIVITY  GRADIENT  PROFILES 

50  X1=X1/X0IV 

yi=yti>/ydiv 

Y2=rtIfl»/YOIV 
IF  fl-l)  52,52,53 

52  CALL  PLOT  (Xi,Yl,3) 

53  CALL  PLOT  (Xi,Yl,2) 

CALL  P-OT  (  XI  «  Y2 , 2) 

54  CONTINJE 

55  CONTINJE 

IF  CNPRO-l)  80,80,55 

SET  UP  A  LOOP  TO  PRINT  AND  PLOT  THE  PIECE-HISE  LINEAR  PROFILES 

56  L=I 0 

60  LP=0 

OO  75  I=1,NDATA 
IF  CHNIIT-AHMAX)  61,61,74 

61  L=LH 
LP=LP*1 
H*S»1 

X (L>-RN(II 
V<LI*HNfI) 

IF  CJPLOT-ll  62,62,55 

PRINT  THE  REFRACTIVITY  PROFILE 

62  X1=X«L) 

Y1=YCL1 

IF  «<RE-)  64,64,63 

63  PRINT  1 4  u  «  L  ,Y-1 ,  XI 

64  IF  INREF)  74,74,70 

PRINT  THE  REFRACTIVITY  GRADIENT  PROFILE 

65  IF  CI-NO ATA )  67,66,56 

66  IF  INPRO-1)  67,67,74 

67  Xi=100  3.MX<M)-XtH,l)  )/(V (H)-Y(N+1)) 
ri  =  fTM)*V{N4U)/2.0 

IF  (KGRAOI  69,69,68 

68  PRINT  140, M, VI, XI 

69  IF  INGRADJ  74,74,70 

PLOT  THE  REFRACTIVITY  AND  REFRA3TI VITY  GRADIENT  PRO-ILES 

70  X1=X1/XDIV 
Y1=YCH1/YDIV 
Y2=YM»l)7YDIV 

IF  ILP-1)  72,72,73 

72  CALL  P-0T  (X1,Y1,3> 

73  CALL  PLOT  (X1,Y1,2> 

CALL  P-OT  (X1,Y2,2) 

74  CONTINJE 

75  CONTINUE 
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C  POSITION  THE  PEN  IF  THERE  HAS  BEEN  A  PLOT 

8C  IF  (JP.OT-1)  90 ,9G,92 
90  IF  (VRE-f  96,96,94 
92  IF  (NGRAD)  96,96,94 
96  CALL  PLOT  (0. 0,0. 0,3) 

CALL  P-3T  (8. 0,-30. J, -31 
CALL  PLOT  (0.C,2.0,-3> 

96  CONTINJE 
RETURN 
C 

110  FORNAT  </// 23 X, * RE F2 ACTIVITY  PROF ILE»/1 OX, *  I* , 8< , 

1*ALTITJ3E*,8X,*REFRACTIVITY»F21X,*(FT)*,11X,*(N-UNITS> *> 

120  FORNAT  ( ///18X , *R£FR ACTIVITY  GRADIENT  PROFILE */l OX ,*I*,8X, 
1*ALTITJQ£*,7X,»REFR  GRADIt NT*/21 X (FT ) * , 9X , * ( N-UNI TS/ K^T) *) 
140  FORNAT  ( 8X , 13, 2(6X , * 1 0 . 2 ) > 

ENO  ' 
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- THIS  ROUTINE  SETS  U3  THE  INTEGRATION  OF  THE  RAY  TRAJECTORY 

- AND  TIME  OF  PR0AGAT10N  EQUATIONS  FOR  EACH  RAY.  WHEN  A  RAY 

- - CROSSES  A  BOUNDARY  3ETHEEN  THE  3IECE-NIS£  LINEAR  SEGMENTS 

- OF  THE  REFRACTIVITY  PROFILE  OR  THE  SOUNOA- t  AT  THE  EARTH'S 

"—-SURFACE.  THE  RAY  EQUATIONS  ARE  INTEGRATED  TO  THE  BOUNDARY 

... - SY  H£AHS  OF  A  VARIA3LE  STEP  SIZE  INTERPOLATION  ALGORITHM. 

- THE  RAY  EQUATIONS  ARE  THEN  RS-IMITIALIZEO  AT  THE  8QUNOARY 

- ftHO  INTEGRATED  TO  THE  NEXT  BOUNDARY.  NHERE  THE  INTERPOLATION 

- --IS  REPEATED. 

COMMON  /ONECOM/  CRH, CRX , CRG.OTR, JAREA , NDATA ,C. Cl, C2 

1  ,RHOC5l,3Q) ,®HI (51,30 I , PI . CMF, CNAUT , AHS 

2  , FRQ, TPM, NHV.NSL.NRMS.ABSRH, PHASE, THETA 

COMMON  /TWOCOM/  RA, DEL X, AHO, ELO ( 51 > , JR A V, EXMAX , AHMA X , XFINAL 
COMMON  /THRCOM/  NEL3, JPLOT , I PLACE , VX ( 3  0 ) .JCASE 

1  ,H151,33) ,G(5l,3G) ,  HN  <50  > ,  RN  ( 50  > ,  S  C  51,  30  > 

2  ,XSIZ,YSIZ,JPLT,XMIN,YHIN,XOIV,YDIV 

3  ,  XM AX , YMAX, NREr, NRA  Y, NPRO, NPLOT.NGRAD 

4  ,<REF,KGRAO,KP_OT 

DIMENSION  YINTtlO) ,3ELXX<5) .PINT tlO) 

DATA  OELXX,NEQ/1000a.,10  0  0.,100.  ,10.,D.  ,2/ 

SET  JP  INITIAL  CONDITIONS 
DO  10  1=1,10 
PINT {11=0.0 

10  YINT (I )  =  S  .  0 
L*1 

ANGLEsO.O 
RHMAS -1 . 0 
STPXsOELX 
VX(LI=CRX 
HCJRAY,L)=CRH 

SIJRAY,L)=ATAN(CRG*RA/CRA*CRH)» 

GCJRAY.LJsPINTdl 
RHO(JRAY,L)=RHMAG 
PHI (JRAY.LI =ANGLE 
YINT  til =CRH 
YINT  C2I =CRG 
VALG=YINTC2> 

VGTEMP=PINT (1) 

VHTEMP=YINT  11 1 
VXTEMP=CRX 
ICOD£s3 
11=0 

CHECN  IF  A  RAY  TRACE  IS  TO  BE  MADE 
IF  INRAY)  14,14,11 

TEST  FDR  MAXIMUM  ALTITUDE 

11  IF  CCRH-AHMAX)  12,12,13 

MAXIMUM  ALTITUDE  NOT  EXCEEDED.  POSITION  PEN  AT  EMITTER  COORDINATES 

12  X=0.0 
Y=CRH>YDIV 

CALL  PLOT  <  X, Y , ICOOE) 
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ICOO* *2 
SO  TO  14 

MAXIMUM  ALTITIOE  EXCEEDED.  POSITION  PEN  AT  UPPER  LEFT  HANO  GRAPH 

13  X=0.0 
T=AHMAX/YOIV 

CALL  P.OT  (X,Y,IC"DE> 

14  IF  (JRAY-1)  15,15,22 

FIND  WHICH  LATER  THE  EMITTER  IS  IN 

15  CONTINUE 

00  20  1=1, NOATA 
IF  (ORH-HN (I)  )  20,21,16 

16  JAREA=I 
JXNTR=I 
GO  TO  25 

20  CONTINUE 
GO  TO  25 
22  JAR£A= JXMTR 
25  CONTINUE 

SET  UP  A  LOOP  FOR  INTEGRATION  Or  THE  ARRAYS  “YINT"  AND  "PINT- 
DO  200  1=2 , IDL ACE  ' 

II=II*1 

CALL  ROUTINE  TO  INTEGRATE  “TINT"  AND  "PINT- 
CALL  R<  (NEQ,CRX,ST»X,YINT,PINT1 

CHECK  HHICH  LAYER  THE  RAY  IS  IN 
IOC  KAR£A=UARE A 

00  120  J=l, NOATA 
IF  (TINT (1 1  -HN ( Jl )  120,120,110 
110  KAREAsJ 
GO  TO  125 
120  CONTINUE 

SET  UP  LAYER  IF  RAY  HAS  INTERSECTEO  EARTH'S  SURFACE 
IF  t YI NT ( 1 > -AHS)  122,122,125 
122  KAREA=N0ATA*1 
125  CONTINUE 

CHECK  WHICH  LAYER  BOUNDARY ,  IF  ANY,  HAS  BEEN  CROSSES  FIRST 
IF  ( JAREA-KARE A J  143,145,130 

AN  UPPER  BOUNDARY  HAS  BEEN  CROSSED 
130  BNSRY=HN(JAREA-11 
KA'REA  =  JAREA-1 
GO  TO  160 

A  LONER  BOUNDARY  HAS  BEEN  CROSSED 
140  BNDRY=HN 1 JARE A) 

karea=jarea*i 

GO  TO  150 

NO  BOUNDARY  HAS  BEEN  CROSSEO,  STORE  ARRAY  VALUES  EVERY 
TENTH  INTEGRATION  STEP 
145  IF  (11*10)  148,146,146 
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146  11=0 
L«LM 
vx( l) =crx 

HCJR4Y,L)=YINT<1> 

S<jRAY,L)=ATAN(YINTC2)*RA/  CRA* YI NT tl> ) > 

GCJR4Y,L)=PINTtl) 

RHO  < JR4Y  ill  =R*MAG 
PHICJR4Y,L)=ANGLE 
148  VALG=YlNTC2) 

VXT£NP=CRX 
VHT£MP=YINT (1) 

VGTEMP=PINTC1) 

ST»X=0ELX 

CHECK  1“  A  RAY  TRACE  IS  TO  BE  MADE 
IF  CNRAY)  200*200,132 

A  RAY  TRACE  IS  TO  BE  MAOE.  CHECK  FOR  MAXIHUM  ALTITUDE  . 

152  IF  CYIMT(l)-AHMAX)  154,154,156 

MAXIMUM  ALTITUDE  NOT  EXCEEDED.  CALL  ROUTINE  TO  PLOT  THE  RAY 
154  X=CRX/XDIV 

Y=YINTC1)/Y0IV 
CALL  °_OT  ( X, Y, ICOOE) 

ICODE=2 
GO  TO  158 

MAXIMUM  ALTITUDE  EXCEEOEO.  TURN  OFF  PLOTTER 
156  ICOOS=3 

158  IF  ( JAREA-KAREA)  173,200,170 

THE  LAYER  BOUNOARV  HAS  BEEN  FOUNO.  SET  UP  FOR  LINEAR  INTERPOLTI OM  SCHEM 
160  OXT3T=l. 

HTEMP=YINT (1) 

CRX=VXT£HP 
YINT  C 1 ) = VHTEMP 
YINT(2)=VALG 
PINT ( 1 )  =  VG  TEMP 

SET  UP  VARIABLE  INTEGRATION  STE3  SIZE  ANO  INTERPOLATE  TO  THE  BOUNDARY 

OO  165  I JK=1»  5 

XX=-VINT(2)/YINTC3) 

XY3HK*<X*XX-2.*(YINT(U-BN0RY»/riNT(3) 

IFCXYCHK.LT.O. )  XYCHK=C. 

YY=SORT (XYCHK) 

CH3X*XK*YY 

IF  {XX.GT.YY1  CHGX=XX-YY 

IF  (I  J<. GT .4)  CHGX=  CBN DRY -YINT C1))7YINTC2) 

SD=$TPX-OXTOT 

IF  CCCH3X.LE.3.) .OR. CCHGX.GE.SD  -  ))  CHGX= 

1  CSTPX-DXTOT) *(BNDRY-YINTC t  > I 7CHTEMP-YINT (III 

CH3X*CHGX-0ELXXCIJK) 

IF  CCH3X.LE.0.)  GO  TO  165 

CALL  ROUTINE  TO  INTEGRATE  TO  THE  BOUNDARY 
CALL  R<  CNEQ,CRX,CHGX, YINT, PINT) 

DXTOTOXTOThCHGX 
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165  CONTINUE 

CHECK  IF  A  RAY  TRACE  IS  TO  BE  MADE 
IF  INRAY-i)  170,152,152 

CHECK  IF  RAY  HAS  INTERSECTED  EARTH'S  SURFACE 
170  IF  <*AR£A- (NDATA+1) >  190,180,183 

RAY  HAS  CROSSED  ZERO  ALTITUDE  BOUNDARY.  FIND  INCIDENT  GRAZING  AMGL 
180  THETA=A3S( ATAN ( YINT (2 > >> 

CALL  ROUTINE  TO  CALCULATE  COMPLEX  SCATTERING  COEFFICIENT 
CALL  SCATT 

SET  UP  COMPLEX  SCATTERING  COEFFICIENT 

RMMAG=RHMA6.*ABSRH 

ANGLE= ANGLE ♦PHASE 

ADD  MULTIPATH  RAY 
OO  185  LL=2 ,13 
185  YINT(LLI=-YINT(LL) 

KAREA  =  NOAT  A 

SET  UP  TO  INTEGRATE  FROM  THE  BOUNDARY  TO  THE  NEXT  “0ELX“ 

190  JAREA=<AREA 
VXTENP=CRX 
WMTEMP=YINT  til 
VAL5- YINT ( 2 ) 

VGTEMP=PINT (1) 

ST«'X*STPX-OXTDT 

CALL  ROUTINE  TO  INTEGRATE  FROM  THE  BOUNDARY  TO  THE  NEXT  “DELX" 

CALL  R<  INEQ, CRX,ST’X, YINT ,PINT) 

CHECK  FOR  MORE  BOUNDARY  CROSSINGS 
GO  TO  103 
200  CONTINJE 

TURN  O"  THE  PLOTTER  IF  THERE  IS  A  RAY  TRACE  PLCT 
ICOOE-3 

IF  CNRAYI  233,230,210 
210  CALL  P.OT  ( X, Y , ICOOE I 

CHECK  IF  THIS  IS  A  PLOT  OF  THE  LAST  RAY.  IF  SO, 

POSITION  THE  PEN  FOR  THE  NEXT  PLOT. 

IF  t JRA V-NELO)  230,220,220 
220  CALL  PLOT  (0.0,0.0,31 

CALL  P.OT  116.0, -30.0, -3) 

CALL  PLOT  (0. 0,2. 0,-31 
230  CONTINUE 
RETURN 
END 
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SUBROUTINE  RK  (N,XN,H,V,P) 

- THIS  ROUTINE  INTEGRATES  THE  RAY  TRAJECTORY  AND  TIME  OF 

- PROPAG'TION  DIFFERENTIAL  EQUATIONS  BY  MEANS  OF  A  FOURTH 

— - — ORDER  RUNGE-KUTTA  A.GORITHM  USING  AN  INTEGRATION  STEP 

- SIZE  0r  H  =  XDELTA/1G ,  WHERE  XDELTA  IS  THE  DISTANCE  INTERVAL 

- - ALONG  THE  EARTH'S  SJRFACE  FOR  PRINTING  AND  PLOTTING  THE 

— - HEIGHT  GAIN  CURVES,  OR  H  *  CHGX  WHERE  CHGX  IS  A  VARIABLE 

- STEP  SIZE  SET  BY  THE  INTERPOLATION  ALGORITHM  IN 

- SUBROUTINE  RKINTG. 

DIMENSION  Y (10),P(1I) ,YODT  (10 > ,PDOT 1 10 > , Q(10 , 4> , R t ID ,4) 
1,YN(10» ,PN(10> 

SET  UP  INITIAL  CONDITIONS 
OO  5  1=1, N 
rNtn=Y<ii 
5  PN(I)=»(IJ 

SET  UP  A  LOOP  TO  INTEGRATE  THE  DIFFERENTIAL  EQUATIONS 
OO  60  1=1,4 

CALL  ROUTINE  TO  SET  UP  THE  DIFFERENTIAL  EQUATIONS 
CALL  ATMOS  (Y,P, YOOT , POOT) 

INTEGRATE  THE  DIFFERENTIAL  EQUATIONS 
GO  TO  <10, 20, 30, 40), L 
10  OO  15  <= 1 , N 

Q(K, L ) *H*YDOT <K) 

15  rm=YN<KI  ♦Q(K,L>/2. 

Rfl,Ll=H*POOT(l) 

P(H=PN(ll+R(l,L)/2« 

XrXN*H/Z. 

GO  TO  5C 

20  OO  25  <*1, N  . 

Q(K«LJ=H*VOOT (K) 

25  Y(<I*YN(KI  +  Q (< , L ) /2. 

Rtl,L)=H*PDOT(ll 

P(ll=PN(l)*R(l,Ll/2» 

X*XN*H^2. 

GO  TO  50 
30  OO  35  <=1,N 

Q(K,LI  *H*YOOT (K» 

35  T (K)=YN(K)+Q(K,L) 

R(1 ,  L) *H*PDOT (II 
P(ll*PN(ll  *R(1,L) 

X=XN+H 
GO  TO  50 
40  DO  45  <=1,N 

Q(<,L) = H* YOOT  <KI 

45  Y(R)*YN(KM(Q(K«1)+2,0*Q(K,2)*2,Q*Q(K«3)  'Q(K,4ll/6.3 
RCl,L»*H*POOT(l) 

.  Ptl>=PN(l)«(R(l,l)+2.0*R(l,2>+2.0»R(l,3>»R(l,4>)/6.0  .. 

XN=  XN*H 
50  CONTINUE 
60  CONTINUE 
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CALL  ROUTINE  TO  FIND  THE  VALUES  OF  THE  DIFFERENTIAL  EQUATIONS 
AT  THE  ENO  OF  THE  INTEGRATION  STEP 
CALL  ATMOS  (Y , P, TOOT , PDOT > 

STORE  THE  NEW  DERIVATIVE  OF  THE  RAY  SLOPE 

V  U)  *  YDOT  (21 

RETURN 

ENO 
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SUBROUTINE  ATMOS  (Y,P, YOOT,PDOT) 

C 

C - THIS  ROUTINE  COMPUTES  ATMOSPHERIC  REFRACTIVITY  AT  THE 

C - RAY  ALTITUDE  AND  SETS  UP  THE  RAT  DIFFERENTIAL  EQUATIONS 

C 

CONDON  /ONECOM/  CRH, CRX, CRG, OTR, JAREA. NO ATA,C, Cl, C2 

1  ,RH0(5l,30),PHI (51,30)  ,  ®I ,  CNF,  5NAUT  ,  AHS 

2  ,*?Q,TPH,NHV,NSL,NRMS,A3SRH, PHASE, THETA 

COMMON  /THOCOM/  RA, DELX, A  HO, ELO< 51 > , JRA Y , EXMAX, AHHAX, XFINAL 
COMMON  /THRCOM/  HELD, JPLOT , 1PL ACE , VX (30 > , JCASE 

1  ,H (51*301 ,G (51,301 ,HN (50) , RN ( 50 > , S 151,  30 ) 

2  *  XSIZ , YSIZ, JPLT ,XMIN,YMIN,XOIV,YOIV 

3  *  XMAX, YMAX,NREr, NRAY, NPRO ,N PLOT, NG RAO 

4  ,<REF,KGRAO,KPlOT 

DIMENSION  Y  (10) ,P (13 ) , YD3T (10 ) , ®00T (10 ) 

C 

C  SET  UP  INITIAL  CONDITIONS 

CM* Y ( 1 ) 

C6=  Y ( 2) 

C 

C  TEST  FDR  THE  APPROPRIATE  ATMOSPHERIC  MODEL 

IF  (NPRD-1)  10,10,23 
C 

C  FREE  SPACE  ANO  EXPONENTIAL  MODELS 

10  REFR=Cl*EXP(-C2*CH> 

OLNDH=-C2*REFR*1.0£-06/(RSFR*i.0E-06  ♦  1.0) 

60  TO  100 
C 

C  THE  ATMOSPHERE  IS  STRATIFIED.  SELECT  THE  APPROPRIATE  MOOEL 
20  IF  (JAREA-1)  SO, 50, SO 
C 

C  ALTITUDE  IS  ABOVE  THE  HIGHEST  REFRACTIVITY  PROFILE  DATA  POINT 
C  USE  AN  EXPONENTIAL  MOOEL  WHICH  "ITS  THE  DATA 

50  REFR*C1»EXP(-C2*CH) 

DLNDH*-C2*REF$*1.0E-06/IREFR*1.0E-0&  ♦  1.0) 

GO  TO  IOC 
C 

C  ALTITUDE  IS  BELOW  THE  HIGHEST  REFRACTIVITY  PROFILE  DATA  POINT 

C  USE  A  »IECE-WIS£  LINEAR  MOOEL 

60  SLOPE* (RN( JAREA-1) -RN( JAREA)  )/HN( JAREA-1) -HN( JAREA)) 
B*RN(JAREA)-SLOPE*HN (JAREA) 

C 

C  COMMUTE  REFRACTIVITY  FOR  THE  PIECE-WISE  LINEAR  NODE. 

SO  REFR=SLOPE*CH  ♦  9 

OLNOH*(SLOPE*1.0E-03)/(REFR*1.0E-06  +  1.0) 

C 

C  COMPUTE  THE  DERIVATIVES  FOR  RAY  TRACES  • 

100  RAD*RA*CH 

YD3T(2)*RAO*(2.0*CRA*CG/RAO)**2  ♦ 

1  DLNOH*RAOM(RA*CG/RAD)**Z  ♦  1,0)**2  ♦  1.0)/R#**2 

TD3T(1)*CG 
C 

C  .  COMPUTE  THE  DERIVATIVE  FOR  TIME  OF  PROPAGATION 
PDDT(1)=RAO*CREFR*1.0E-06*1.0) 

1*SDRT (1 . 0 ♦ ( RA*CG/RAO>  **2) / (C*RA) 

RETURN 

ENO 


Fig.  6l.  (continued) 
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SUBROUTINE  SCATT 

- THIS  ROUTINE  CALCULATES  THE  COMPLEX  SPECULAR  SCATTERING 

- COE'FIOIENT  FDR  WAVEFRONT  REFLECTION  FROH  SMOOTH  OR  ROUSH 

- LANO  AND  SEA  SURFACES 

COMMON  /ONECOM/  CRH, CRX , CRG.DTR, JAREA, NOATA ,C, Cl , C2 

1  , RH0<51, 30),PHI(51,30),PI,CMF,CNAUT,AHS 

2  ,-R0,TPW,NHV,NSL,NRMS,A9SR-t,  PHASE,  THETA 
DIMENSION  £PS(3,2),3IGMA(3,2),OELHSL(9,2) 

REAL  MI, NR 

DATA  E»S,IC/8C .,69. ,65., 4. ,10 .,30.,D/ 

DATA  SI3HA/L.3,6.5,lo.O, . 0G01, .3 016,.  01/ 

DATA  DELHSL/Q  « , ,2, , 6, 1,1 « 1 ,7, 2. 5, 4,3*8. 6, £2,9 
1 1  0. «9«  «  3G«  < 56* 1 112*  «  214« »429« ti268*  « 2146*/ 

• 

SET  UP  INITIAL  CONDITIONS 
IF  (ICI  10,10,50 
10  IC*1 

JJ=MRHS+1 

IF  (NSL-ll  12,12,30 
12  KK=1 

IF  (FRD-lSOOOCOOOO..)  14,16,16 
14  11*1 

GO  TO  40 

16  IF  (FR0-50 OOOOOOOO,)  18,20,20 
18  11*2 

GO  TO  40 
20  11*3 

GO  TO  43 
30  KK*2 

IIsNSL-1 

CALCULATE  THE  COMPLEX  DIELECTRIC  CONSTANT 
40  ER=EPS(  1 1 ,  KK) 

EI=60.0*SIGMA(II,KK)*C/(CMF*FRQ) 

a=sort (er**2+ei**2) 

ALPHA*ATAN2  <£I,ER) 

OELTAH*OELHSL ( JJ.KK) 

CHECK  WHICH  POLARIZATION  IS  BEING  USED 
50  IF  (NHtf-1)  60,60,70 

HORIZONTAL  POLARIZATION 
60  NR=SIN ( THETA) -S5RT (A)* CO S( ALPHA/2,) 

NI*SQRT (A) *SIN (ALPHA/2.) 

DR=SIN( THET  A) *S0RT ( A I *COS (ALPHA/ 2 , ) 

OI*SQRT ( A) *SIN (ALPHA/2,) 

GR*(MR*DR-NI*0I>/(0R**2OI**2> 

GI*  f NR*D1*NZ*0R) /(0R**2*0I**2) 

GO  TO  SO 

VERTICAL  POLARIZATION 

70  NR*SaRr(A)*COS(ALPHA/2.)*SIN(THETA)-i.O 
NI*S0RT(A)*SlNCAL°HA/2.)  *SIN(THETA) 

OR*SQRT (A) * COS (ALPHA /2. )*SIN( THETA) *1,0 
DI*SDRT (A) *SIN( ALPHA/2. )*SIN( THETA) 


Fig.  6l.  (continued) 
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GR*<MR*DR+NI*0I»/<0?**2*0I**2» 

GI*«NR*0I-NI*DR)/fD?**2«-DI**2> 

CALCULATE  THE  COMPLEX  FRESNEL  REFLECTION  COEFFICIENT 
60  A3SRH*SaRT<GR**2*GI*»2l 
PHAS£*ATAN2(>GZ*GRI 

CALCULATE  THE  TOTAL  SPECULAR  REFLECTION  COEFFICIENT 
ABSRHsABSRH*EXPt-0.5*  <«».  0*PI*DELTAH*SIN  ITHETA) •FRQ/C>**2) 
RETURN 
ENO 


Fig.  61 .  (continued) 
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SUBROUTINE  HGAIN 

- THIS  ROUTINE  COMPUTES  RELATIVE  rIELD  STRENGTH  AMO  POWER 

- DENSITY  (NORMALIZED  TO  FREE  SPACE  VALUES).  HEI3HT  SAINS 

- - ARE  OBTAINED  FROM  CALCULATIONS  AT  100  “WINOOWS“  OR  ALTITUDE 

— - INCREMENTS  EXTENOIN3  VERTICALLY  BETWEEN  THE  HIGHEST  AND 

- --LOWEST  RAYS. 

COMMON  /ONECOM/  CRH, CRX , CRG , DTR, J AREA , NO AT A , C , Cl , C2 

1  ,  RHD(5l,30),PHI(51,30> , PI , CMF, CNA UT, AHS 

2  ,CRQ,TPW,NHV,N5L,NRMS,ASSRH, PHASE, THETA 

COMMON  /TWOCOM/  RA, OELX, AHO,£LO( 51) ,JRAY, EXMAX, AHMAX,XFINAL 
COMMON  /THRCOM/  NELO, JPLOT , IPLACE , VX ( 3 0 > , JCASE 

1  ,M(51«3C),G(51,3u)«HN(5D),RN(53),S(51,30) 

2  ,  XSIZ,YSIZ,JPLT»XMIN, YMIN, XDIV, YDI V 

.  3  i  XMAX,.YHAX,NREF,NRA  Y  ,NPRO,NPLOT ,  NGRAD 

4  , KRtF , KGRAO, KP.OT 

DIMENSION  L (51),T0IFF(3ai ,ESIG(30 I ,PSC(3C ) 

1,T3IG(3C) , ESIGT(30).PSCT(3C),TSI6T(30) 

DATA  ZDIV/40.0/ 

SET  UP  INITIAL  CONDITIONS 

Rl*RA*AHO 

XN-DELX/XOIV 

YN* AHMAX/YDIV  . 

VS* AHSfYOIV 
NP*100 

SET  UP  A  LOOP  TO  CALCULATE  RELATIVE  FIELD  STREN5TH  OR 

POWER  DENSITY  AT  EACH  INCREMENT  OF  OISTANCE 

DO  1Q0G  K=2, IPLACE 

HL0*AHMAX 

U=VX ( X) /RA 

FIND  HIGHEST  AND  LOWEST  RAYS 
DO  IOC  1=1, NELO 
HLO*AMINl (HLO, H ( I , K) ) 

IF  (HLO.EQ.H(I,KI)  IRA*I 
100  CONTINUE 

HHI*AHI.M1(H  (1,K),AHMAX) 

MDIFF*WHI-HLO 
HTB*A83 (H(1,K>  -H  (2,0  ) 

EXIT  IF  ALL  RAYS  EQUAL  OR  EXCEED  THE  MAXIMUM  ALTITUDE 
IF  IHLO-AHMAX)  200,1100,1100 

SET  UP  WINDOW  SIZE  AND  WINDOW  POSITIONS  TO  EXCLUDE  THE 
HIGHEST  ANO  LOWEST  RAYS 
200  WINDOW* HOI FF/95. 

MNOW*HMI -WINOOW/IO , 

hfimal*hlo*minoow/io. 

RWP*HNOW-HFINAL 

FIND  THE  ALTITUDE  INCREMENT  FOR  EACH  NEW  WINDOW  POSITION 
DH*RWP/FLOAT(NP) 

HNOW*HMOW*DH 


Fig.  6l  .  (continued) 
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C  SET  UP  A  L00O  TO  POSITION  THE  WINDOW  EVERT  •*DH"  FEET  IN  ALTITUD 
03  900  J=1,NP 
HNOW=HNOW-OH 

CALCULATE  THE  UPPER  ANGLE  LIMIT  SUBTENDED  BY  THE  WINDOW  IN 
FREE  S»AC£ 

R2  =  RA-HNOW 

DSQsRl**2  ♦  R2**2  -  2. 0*R1*R2*C3S ( U> 

E1=ACD3 ( (OSQ  ♦  <R1*TAN(U) ) »»2  -  ( ff 1/COS ( U)-R2> ** 2> 
1/<2.0*R1»TAN(U)*SQRT (OSQ) ) ) 

IF  CR1/COSCU) .GT.R2)  El=-£1 

CALCULATE  THE  LOWER  ANGLE  LIMIT  SUBTENOED  BY  THE  WINDOW  IN 
FREE  SPACE 
R2=RA«-HN0W-WIN00W 

DSQ=R1**2  +.  R2**2  -  2.0*R1*R2*C3SIU) 

E?= ACOS ( (OSQ  ♦  (Rl’TANIU) » **2  -  < R1/C0S ( U) -R2 ) **2> 

1/(2. 0*R1*TAN(U)*S(3RT< OSQ))  ) 

IF  (R1/C0S(U) .GT.R2)  E2*-£2 
EDIFr =SQRT ( ABS (E1-E2) /DTR) 

CALCULATE  THE  ELEVATION  ANGLE  O'  THE  EMITTER  AT  THE  WINDOW 
Rl  =  RAHNOW-WINDOW/2. 

R2=RA* AHO 

DSQ=R1**2  ♦  R2**2  -  2.0*R1*R2*C3S(U> 

AO=ACOS((OSQ  ♦  (R1*TAN(U) > **2  -  ( Rl/COS CU> -R2) **2> 
1/(2.0*R1*TAN(U)*SQRT(DSQ) ) ) 

IF  (Rl/COS(U) .GT.R21  A0=-A0 

SET  UP  THE  WINDOW  ALTITUDE  PLUS  ITS  UPPER  AND  L-WER  BOUNDARIES 
HU»sHNOW 

HH=HNOW-WINOOW/2« 

HDN=MNOW-NINDOW 

SET  UP  CONTROL  INTEGERS  FOR  EACH  RAY 
DO  220  I  =  1 . NE L 0 
IF  (H(I, K) .GT.HUP)  SO  TO  210 
IF  (H(I,K) .LT.HON)  SO  TO  215 
L«)=2 
GO  TO  220 
210  LCIl =3 

GO  TO  220 
215  L(I»=1 
220  CONTINJE 

SET  UP  INITIAL  CONDITIONS 
E=0.  Q 
JSLOPE=C 
JREV=*0 
USTART  =  0 
JSTOP»0 
LC  =  0 

TNAXsC.O 
KT3T=Q 
MEL*NE.3-1 
OO  222  IFS=1»  2 
TDIFF (IrS)=0.0 


Fig.  6l ,  (continued) 
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222  ESISd^SlsO.O 

IF  (L(l» .NE.2I  GO  T3  225 

INITIAL  CONDITIONS  IF  FIRST  RAY  IS  INSIDE  THE  WINDOW 
JST  ART  =  1 
LC=1 
Jl  =  l 

225  CONTINJE 

SET  UP  A  LOOP  WHICH  SCANS  RAY  HEIGHT  VERSUS  ELEVATION  ANGLE 
00  600  1=1, HEL 
230  JLAST  =  JSL0PE 

JC=IAB5(L(I)-L(I+1>) 

IF  IHI)-L(MI)  260,400,  240 

NORMAL  RAY  DRDER 
240  JSLOPE=l 

IF  (JLAST . EQ.2)  GO  TO  310 
GO  TO  270 

INVERTED  RAY  ORDER 
260  JSLOPE*? 

IF  (JLAST.EQ.il  GO' TO  310 

CHECK  IF  NEXT  RAY  IS  THE  LOWEST 
270  IF  (I*1>IRA)  200,272,280 

NEXT  RAY  IS  THE  LOWEST.  CHECK  IF  IT  HAS  BEEN  REFLECTED  FROM  THE  EARTH 
272  IF  (RHD (IRA *K) -1.0)  274,290.290 

LOWEST  RAY  HAS  BEEN  REFLECTED.  SET  ITS  IMAGE  ALTITUOE  BELOW  THE  EARTH 
274  H(I»1,<I=AHS-H(X«1,<> 

GO  TO  290 

CHECK  IF  THIS  RAY  IS  THE  LOWEST 
200  IF  (I-IRA1  290,600,290 

EITHER  ONE  OR  TWO  WINOOW  LIMIT  CROSSINGS 
290  IF  (JC-l)  360,300,300 

ONE  WINDOW  LIMIT  CROSSING 
300  LC=LC+1 

IF  (LC-l)  340,340,350 

RAY  ORDER  REVERAL  OCCURS.  CHECK  IF  RAY  BUNOLE  EXISTS  INSIOE  WINDOW 
310  IF  (JC-ll  320,320,230 

CHECK  IF  RAY  REVERSAL  OCCURS  INSIDE  OR  OUTSIOE  THE  WINOOW 
320  IF  (LC)  230,230,330 

RAY  ORDER  REVERSAL  OCCURS  WITHIN  THE  WINDOW 
330  J2= I 
JREV*1 
GO  TO  440 

FIRST  WINOOW  LIMIT  CROSSING 
340  J1=I«1 


Fig.  6l .  (continued) 
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J2=I*1 

IF  (X.EQ.MEL.AND.L(NELO)  . EQ. 2)  30  TO  420 
GO  TO  300 

SECOND  WINDOW  CROSSING 
360  J2=I 

GO  TO  ‘*40 

BOTH  HINOOH  LIMITS  CROSSED .  CHECK  FOR  NEARBY  C4JSTICS 
38C  IF  CABS<H<I«K)-HmitKn-5.0*HT3>  390,390,560 

BOTH  HINOOH  LIMITS  CROSSED.  NO  CAUSTIC  HAS  BEEN  FOUNO 
390  J1=I*1 
J?=I 

GO  TO  440 

NO  HINOOH  LIMIT  CROSSINGS 
400  IF  CI.Ea.NEL.ANO.L(NELO) .EQ.21  GO  TO  420 
GO  TO  SOC 

LAST  RAT  IS  INSIDE  WINDOW 
420  J2=NEL3 
JST3P=1 

THE  INCREMENT  OF  RAYS  AT  THE  WINDOW  HAS  BEEN  FOUND. 

CHECK  IF  RAY  ORDER  IS  NORMAL  OR  INVERTED. 

440  IF  (JLAST-l)  460,463,469 

NORMAL  RAY  ORDER.  CALCULATE  THE  UPPER  AND  LOHER  ANGLES 
SUBTENDED  BY  THE  HINOOH. 

460  EU=»  =  EL3(J1> 

EDN=ELOIJ2» 

IF  USTART.EQ.O)  EU** 

lEL3<Jl-U-<ELO<Jl-ll-EL3<Jl>>*tH{Jl-i,K>-HUP)/<-l(Jl-l,IO-H(Jl,IO> 
IF  IJST3P.EQ.0.ANO. JREV.EQ.O)  £DN  = 
lELOU2J-(ELO(J2»-EL3(J2+l))MH(U2,IO-HDN>f(HC  J2, K) -H C J2*i, Kl > 
E=S3RT (EUP-EON) 

PSCAT=»HI(J2,KJ 
RSC  AT  =  RHO  ( J2,  Kl 
GO  TO  500 

INVERTED  RAY  ORDER.  CALCULATE  THE  UP»ER  AND  LOHER  ANGLES 
SUBTENDED  BY  THE  HINOOH. 

480  E0N=SL3CJ1> 

EUPsELD  C J2 ) 

IF  (JSTART.EQ.C)  E0N= 

1ELD<J1-1»-  (EL0(J1-1>  -ELO(  Jl)  J*H(  J1 -1,  K) -HDN) /  HI  J1  -1,K)-HCJ1,  K> » 
IF  UST3P.EQ.O.ANO.JREV.EQ.OI  EUP  = 

IEL3<J2)-»EL0<J2J-EL3<J2*1>  »  MHU2  ,  K» -HUP>  H ( J2,  K> -H CJ2+ 1, Kl  » 
E-S3RT  C  EDN-EUP) 

PSSATs»HI( J1,K) 

RSCAT=RH0(J1*K» 

SET  UP  INITIAL  CONDITIONS  TO  CALCULATE  MEAN  ELEVATION  ANGLE 
AND  MEAN  TIHE  OF  ARRIVAL 
500  KTOT*KTOT+i 
JA  =  J1 


Fig.  61  (continued) 
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JB=  J2 
JCOUNT®0 

IF  CJ1-J2)  528.520.310 
510  JA  =  J2 
JB=  J1 

CALCULATE  MEAN  ElF>'ATI0N  ANGLE  AND  TIME  OF  ARRIVAL  FOR  EACH 
WAVEFRONT 
520  SN£W=0.0 
TNEW=0.C 

00  540  JAB=JA, JB 
JCOUNT  =  JCOUNT  *-1 
SNEW=SNEW*S  (JAB. K) 

540  TNEW=TNEW*G  (JAB.K) 

ALPH=-SNEW/FLOAT (JCOUNT) 

TSIGTf<TOT)=TNEW/FL0AT (JCOUNT) 

CALCULATE  ANGLE  INCREMENT  OF  EACH  WAVEFRONT 
ESI GT (<TOT)=E*RSCAT*COS(AO )/COS(ALPH) 

PSCT (KTOT) =PSCAT 

RESET  INITIAL  CONDITIONS  TO  CONTINUE  SCAN  OF  HEIGHT  VS  ANGLE 
560  H(IRA.O=ABS(H(IRA,<)} 

LC  =  0 
JREV=0 
JST  ART  =  0 
600  CONTINJE 

CHECK  IF  NINOOH  IS  IN  A  SHADOW  REGION  RESULTING  FROM  A  CAUSTIC 
IF  (KTOT)  780,780,610 

WINDOW  IS  NOT  IN  A  SHAOOW  REGION.  *>UT  WAVEFRONTS  IN  ORDER 
0*  THEIR  TIMSS-OF-ARRIVAL 
610  CONTINJE 

DO  660  I S® 1  .KTOT 

TMINslj.O 

DO  650  IR=1 .KTOT 

THIN* AMIN1 (THIN, TSI 3T (IRI ) 

IF  (TMIN.EQ.TSIGTdR))  ININ=IR 
650  CONTINJE 

TSIG(IS) =TSIGT (I MIN) 

ESIG(IS)=ESIGT(IM1N) 

PSC(IS)=PSCT(IMIN) 

TSIGT(ININ)=20.0 
660  CONTINJE 

CHECK  IF  MORE  THAN  ONE  WAVEFRONT  IS  PRESENT  IN  THE  WINDOW 
IF  OCTOT.EQ.l)  GO  TO  700 

CALCULATE  TIME  DIFFERENCE  OF  ARRIVAL  (T.D.O.A.)  BETWEEN  WAVEFRONT 

DO  670  IT =2,KT0T 

TOIFFtlT)sTSIGCIT)-TSIG(l) 

670  THAKsANAKKTMAX.TOIFFdTI) 

SET  UP  A  LOOP  TO  CALCULATE  TOTAL  FIELD  STRENGTH  OR  »OWER 
OENSITT  IF  WAVEFRONTS  OVERLAP  IN  TIME 
700  ESUNRsO • 0 


Fig. 


61.  (continued; 


non  o o  o  o  no  oo  oo  oo  r> o  n  o 


162 


ESUHI=0.C 

DO  720  IFI=1, KTOT 

IF  (TDIFF(IFI» ,GT.T»HJ  GO  TO  71J 

ES'JMR=ESUMR«-ESIG(IFn*COS  CPSC  (IFI)+2.G»PI*FRQ*TDIFF(IFI)  ) 
ESUNI=ESUMI*£SIG (IFI ) *SIN(PSC(I”I>*2 . 0*PI*FRQ*T0 IFF(IFI)) 

710  CONTINUE 
720  CONTINUE 

£SUN=S3RT  (£SUNR**2+2SUMI**Zi 

CHECK  IF  FIELD  STRENGTH  OR  POWER  DENSITY  IS  TO  BE  ODHPUTED 
IF  CNPLOT-1)  730,743,760 
730  IF  (KPLOT-l)  740,740,760 

CALCULATE  RELATIVE  FIELD  STRENGTH  OF  EACH  WAVEFRONT 
740  CONTINUE 

DO  750  IFS=1, KTOT 

750  ESIG  CI::S)  -10.  *ALOG10  (ESIGdFSI/EDIFF) 

CALCULATE  TOTAL  RELATIVE  FIELD  STRENGTH 
SP*10 ,*ALOGlQ  (ESUH/EOIFF) 

GO  TO  900 

calculate  relative  power  oensity  of  each  wavefront 

760  CONTINUE 

DO  770  IFS=l,KTOT 

770  ESIGCI-SI=20.*ALOG13 (ESIG ( IFS1 /•  OIFF) 

CALCULATE  TOTAL  RELATIVE  POWER  DENSITY 
SP=2D.*ALOG10 (ESUH/EDIFFI 
GO  TO  500 

SET  FIELD  STRENGTH  AND  POWER  OENSITY  FOR  A  SHADOW  REGION 
760  SP=-103GD0.0 

ESI G(l) -'100000,0 
ESIG(2)*-1000U)0.0 

PRINT  HEADING  AND  RAY  HEIGHTS  I-  THIS  IS  A  NEW  DISTANCE 
600  IF  (J.NE.l)  GO  TO  820 
XNAUT*VX(K)/CNAUT 
IF  (KPLOT-l)  616,812,814 
612  PRINT  1200 , XNAUT 
GO  TO  916 

814  PRINT  1210, XNAUT 

POSITION  THE  PEN  AN0  DRAW  THE  NEW  ORDINATE  AXES  IF  THERE  IS  A  P.OT 
816  IF  (NP'.OTJ  820,620,918 
818  CALL  PLOT  <XN,0.Q,-3) 

CALL  PLOT  (0.0, YS, 3) 

CALL  P.OT  (0. 0 ,YN,2) 

CALL  P.OT  (0.0, 0.0, 3) 

820  CONTINUE 

.  PRINT  THE  RELATIVE  FIELD  STRENGTH  OR  POWER  DENSITY  ANO 
TINE  DIFFERENCES  OF  ARRIVAL 
IF  (KPLOT)  826,826,922 

822  PRINT  1220,  J,  HH,  SP,  ESIGdl  ,  ESIG(  2  )  ,TMA  X  ,  TDI FF  ( 2)  ,  KTOT 
IF  (KTOT ”1 )  826,826,824 
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824  PRINT  1230, II ESIG ( IFS) ,TOIFF{ IFS J > ,IFS=3,KTOT> 

C 

C  SET  PARAMETERS  IF  THERE  IS  A  PL3T 
826  IF  INPLOT)  900.900,350 
650  Y=HH/YDI V 

IF  IJ-l)  854,854,855 
C 

C  SET  FIRST  END  POINT  OF  PLOT 

854  X-0.0 
ICODE-3 

CALL  P.OT  ( X, Y, ICOOE) 

C 

C  CHECK  IF  FIELD  STRENGTH  OR  POWER  DENSITY  IS  TOO  LOW  FOR  PLOTTING 

855  IF  ISP* 15. 0  )  856,857 ,857 
C 

C  FIELD  STRENGTH  OR  POWER  OENSITY  LESS  THAN  -40  03.  TURN  OFF  PLOTTER. 

856  XS-15.0/ZDIV 

CALL  PLOT  ( X, Y, ICODE ) 

ISODE=3 
GO  TO  358 
C 

C  PLOT  TOTAL  RELATIVE  FIELD  STRENGTH  OR  POWER  DENSITY 

857  X=S°/2DI V 
ICOOEs? 

CALL  PLOT  ( X, Y .ICOOE ) 

C 

C  SET  END  POINT  OF  PLOT 

858  IF  I J- M° )  900,859,859 

859  X=0.0 
IC30E=2 

IF  ISP. LT, -15. 01  IC00E=3 
CALL  PLOT  IX, Y, ICODE) 

900  CONTINUE 
1000  CONTINUE 
C 

C  POSITION  THE  PEN  IF  THERE  HAS  BEEN  A  “LOT 

1100  IF  INPLOT)  1120,1120,1110 
1110  CALL  P.OT  10.0,0.0,31 

CALL  PLOT  18. 0,-30. 3, -3) 

CALL  PLOT  10.0,2.0,-31 
1120  CONTINUE 
RETURN 
C 

1200  FORMAT  l /✓/, IX, 'PRINTOUT  OF  FIELD  STRENGTH  AND  T.O.D.A.'S  AT  X  =' 
l, F8. 2, 2X, 'NAUTICAL  NILES', 

2/F3X,*ND.*,4X, 'HEIGHT  |FT)*,4X,*  TOTAL  FIELD  (03)  ',4X, 

3*OIRECT  FIELD  tOB)  *, 4 X,* MULTIPATH  FLO  103)', 4X, 

4'HINIMJM  TOOA  tSEC ) » , 4X , 'TDOA  13EC)*,4X, 'NS'/) 

1210  FORMAT  (///, IX , 'PRINTOUT  OF  POWER  DENSITY  AND  T.D.O.A.'S  AT  X  =* 

1, F6. 2, 2X, 'NAUTICAL  MILES', 

2F/3X, 'NO. *,4X, 'HEIGHT  (FT) ',2X, 'TOTAL  POWER  DENS  ID3)',2X, 

3'OIREDT  POWER  OENS  (OB) ', 2X, 'MULTIPATH  PWR  (D3)',4X, 

.  4*MINIMJM  TOOA  (SEC )', 4X, 'TOOA  ( SEC) ', 4X , *NS*/ ) 

1220  FORMAT  ( 2X, 13, 5X.F9. 2, 9X , 3 (F9 .4, 1 3X) , 1PE11 . 4, 1PE18 . 4,4X, 12) 

1230  FORMAT  I 72X ,F9. 4, 24X , 1PE18 . 4) 

ENO 


?i£.  61.  (continued) 


SUBROUTINE  PLOTTR 
C 

C - THIS  ROUTINE  SETS  UP  THE  AXES  FOR  ALL  PLOTS 

C 

COMMON  /TWOCOM/  RA , OELX , AMO, ELO ( 5 1 ) , JR A Y , EXMAX , AHHAX .XFINAL 
COMMON  /THRCOM/  NELO. JPLOT , IPLACE , V X ( 30 > , JC ASS 

1  ,-1(51,33)  ,G<51,30>  ,HN(50)  ,RN(50>  ,S  151,30) 

2  ,XSIZ,TSIZ, JPLT,XHIN, YMIN,  XOIV.YOIV 

3  , XMAX  ,  YMAX, NREF ,  NRAY,NPRO,MPLOT, NGRA'D 

4  ,<REF,KGRAD»KP_OT 

DIMENSION  A  HD (8) ,TTL (8 ) , XTL ( 8 ) , Y TL f 8)  ,  TLE ( 8 )  ,  ZT. tl) 

DATA  YNIN, YSIZ, SIZE/0 . 0 ,  5.0,0.150  5/ 

OATA  ZTL ( 1 ) *IZTL»ZSIZ«ZMIN,ZDIV/1QH  GAIN  (DS),13,1.3 
1,-20.0,40.0/ 

C 

C  SET  UP  INITIAL  CONDITIONS  IF  THIS  IS  THE  FIRST  aLOT 
IF  (JPLT)  10,10,20 
10  YMAX= AHMAX 

YOIV=(YNAX-YNIN)/YSIZ 
TDIVA=YOIV/1QOO.O 
Y3T  sVSI Z  +1 . 0 
Z°T »YSIZ«-0 . 5 
YN= YM A  X/ YOI V 
YS=AHS/YOI V 
JPLT=1 
C 

C  READ  P.OT  SCALE  FACTOR  AND  PLOT  BANNER 
READ  110, SCALE 
REAO  120,IAHD, (AHD(«) ,K=1,7) 

C 

C  SCALE  THE  PLOT  SIZE 

CALL  FACTOR  (SCALE) 

C 

C  POSITION  THE  PEN  AND  PRINT  THE  PLOT  BANNER 
IF  (JCASE-1)  35,15,20 
15  YAHD=(10 .0-SIZE*IAHD)/2. 

CALL  P.OT  (4. 0,-30. 0,-3) 

CALL  P.OT  (0.0, 2. 0,-3) 

CALL  SYMBOL  ( 0 . 0 , YAHO, SIZE , AHO, 30 . 0 , IA HD ) 

C 

C  POSITION  THE  PEN  FOR  THE  NEXT  PLOT 

CALL  P.OT  (6. 0,-30. 0,-3) 

CALL  PLOT  (0.0, 2. 0,-3) 

C 

C  READ  AXIS  LABELS  AND  TWO  LINES  OF  PLOT  TITLES 
20  READ  120.IXTL,  (XTL ( <) » K=l, 7) 

REAO  120.IYTL, (YTL(<) ,K=1,7) 

READ  120,ITTL,  (TTL ( < J , K  =  1 , 7) 

READ  120 , I TLE, (TLE ( <) »K=i,  7) 

C 

C  CHEC<  XHICH  PLOT  THIS  IS 

IF  ( JP.OT-2 )  30,40,50 
C 

C  SET  UP  INITIAL  CONDITIONS  FOR  REFRACTI VITY  PROFILE 
C  AND  GRADIENT  PLOTS 

30  XMAX=4 3  0.0 
XMI N=  0 . 0 


Fig.  6l.  (continued) 


oo  no  non  no  nn  nn  ooo 


60  TO  50 
40  XMAX*100.0 
XHIN*-30C. 0 
SO  XSIZs4 . 0 

XOItf=<XMAX-XMIN)/XSIZ 
XDI VA=X0IV 
60  TO  70 

SET  OP  INITIAL  CONDITIONS  FOR  RAY  TRACE,  FIELD  STRENGTH, 

AND  POWER  DENSITY  PwOTS 
60  XMAX=EXMAX 
XMIN=0 • 0 
XSIZ*10.0 

XDIV=(XNAX-XMIN)/XSIZ 

XOIVA=XFINAL/XSIZ 

XRE=X$IZ-C.S 

SET  UP  PARAMETERS  FDR  AXIS  LABELS  AND  PLOT  TITLES 
70  XTTL-(X5IZ-SIZE*ITTL)/2. 

XTLE=<X3IZ-SIZE*ITLEI/2. 

DRAW  AND  LABEL  THE  AXES 

CALL  AXIS  (0.0,0.0',YTL,IYTL,YSIZ,90.0,YMIN,YDIVA) 

CALL  AXIS  (0.0,0.0,XTL,-1XTL,XSIZ,0.0,XHIN,XDIV4> 

PRINT  TNE  PLOT  TITLES 

74  CALL  SYMBOL  ( XTTL , Y=T,SIZE ,TTL ,3 . 0 , ITTL ) 

CALL  SYMBOL  (XTLE.ZPT, SIZE, TLE, 0.0, ITLE) 

DRAM  THE  REFERENCE  SCALE  FOR  ALL  FIELD  STRENGTH  AND 
POWER  DENSITY  PLOTS 
IF  <J“_3T-31  92,82,30 

00  CALL  AXIS  (XRE,ZPT,ZTL,-IZTL,ZSIZ,G.O,ZHIN,ZOIVI 
CALL  PLOT  f XSIZ, ZPT ,3) 

CALL  PLOT  (XSIZ»YPT,21 
62  CALL  PLOT  (C. 0,0. 0,31 

DRAW  THE  EARTH  SURFACE  IF  NOT  AT  ZERO  MEAN  SEA  LEVEL 
IF  ( AHS)  86,85,64 
04  CALL  PLOT  (0.0,YS,3I 
CALL  PLOT  (XSIZ.YS, 2) 

CALL  PLOT  (0.0,0.0,31 

POSITION  THE  PEN  FOR  THE  REFRACTI VITY  GRADIENT  PROFILE  PLOT 
66  IF  (JPLOT-2)  100,90,100 
90  XNs-XNIN/XOIV 

CALL  PLOT  (XN.O. 0,-31 
100 . CONTINJE 
RETURN 

110  FORMAT  (3F15.7) 

120  FORMAT  (I2,8X,7A10) 

END 
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TA3LE  10 


CARD  INPUT  DATA 


.  CARD. _ VARIABLE _ FORMAT _ OESCRIETION _ _ 


NC4SE  15  NUMBER  OF  PROGRAM  RUMS. 


CODE.  FOR. REFRACTIVITY  MODEL... _ 

NRE F=Q  fREE  SPACE  REFRACTIVITY 

_ HO  DEI _ 

NREF=1  EXPONENTIAL  REFRACTIVITY 

_ MODE! _ 

NREF-2  INPUT  REFRACTIVITY 
_ ERJ3FILE_M00EL _ _ _ 


3  NDATA  15  NUMBER  OF  REFRACTIVITY  PROFILE 

_ DAT  A_LEVSLS_. _ 

NDATA  =  1  IF  NREF=C  OR  1 


RN(L)_ 


_HN  ( 1J _ I S  THE  HIGHEST,  ALTITUDE _ 

IN  THE  REFRACTIVITY  PROFILE. 

RN  (1)  IS  REFRACTIVITY  AT  HN  ( l )_•  __ 
**  (LEAVE  THIS  CARD  BLANK  IF 
_ EITHER  NREFsp  OR  NREF=l» _ 


2F15.7  HN(2>  IS  SECOND  HIGHEST  ALTITUDE 

_  IN  THE  REFRACTIVITY  PROFILE. _ 

RN (21  IS  REFRACTIVITY  AT  HN(2). 

_ **__  (OMIT  CARDS  5  THROUGH  (NDATA±3I 

IF  EITHER  NREF=C  OR  NREF=1» 


NDATA  HN(NOATAl _ 2F15.7  LOWEST  _ALTITUDE  IN  PROFILE. 

♦3  RN(NOATA)  REFRACTIVITY  AT  HN(NDATA) , 

_ '  ♦♦  (OMIT  IF  NREFaQ  OR  NREF=IJL 


NOATA 
♦4 _ 


AHS 

AHD _ 

AHNAX 


3F15.7  EARTH  SURFACE  ALTITUDE  IN  FEET. 

_  EMITTER  ALTITUDE  IN  ' EET .  _ 

MAXIMUM  ALTITUDE  IN  FEET  FOR 
_ _ PRINTEQLAND_PLOTTED  output. 


NDATA  XDELTA  2F15.7  DISTANCE  INTERVAL  IN  NAUTICAL  MI 

♦5 _ _ FOR  PRINT  AND  PLOT  OUT»UT.  _ 

XFINAL  MAXIMUM  DISTANCE  IN  MAUT  MI  FOR 

_ PRINTED  AND  PLOTTED  OUTPUT. 


NDATA  ELOS1  2F15.7  HIGHEST  RAY  ELEVATION  ANGLE  IN 

♦6 _ _  DEGREES. _ _  _  _ 

'LOWEST  RAY_ELEVATION 'ANGLE  IN 
OEGREES. 


EL3S2 
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TABLE  10  (continued) 


(■ 


_**  THE  QUANTITY  (EL0S1-EL0S2) _ 

SHOULD  NOT  EXCEEO  l.C  DEGREE. 


NOATA 
♦  7 


_E M I TTER_ FREQUENCY  IN  MEGAHERTZ . _ 

EMITTER  PULSE  WIDTH  I N  MICROSEC. 
_  SET  PW=1G00300 .0  Is  EMITTER  IS 
CONTINUOUS  WAVE  (C.M.) 


4-8 

NH</  =  1 
NHV=2 

NSL 

CODE  FOR 

NSL  =  1 

NSL  =  2 
NSL  =  3 

NSL  =  4 

NRMS 

CODE  FOR 

HORIZONTAL  POLARIZATION 
VERTICAL  POLARIZATION^ 
EARTH  SURFACE  TYPE. 

SEA  WATER  _ 

VERY  DRY  L"AN0 
"AVERAGE-  LAND 


TABLE  BELOH. 


_ TABLE  OF  EARTH  ROUGHNESS  CODES _ 

(FROM  MAURICE  LONG.  "RADAR  REFLECTION  FROM  LAND  AND  SEA". _ 

REF.  (35),  ANO 

THE  OFFICE  OF  TELECOMMUNICATIONS,  U.S.  DEPT.  OF_CONMER3Ej _ 

REF.  (39)  ) 


CODE  STANOARO  DEVIATIONS  OF  HEIGHT 

_ _ _ _  -OR  SEA_AND  LAND  SURFACES _ • _ 

NRMS  (SEA)  (LANO) 


NDATA 
♦9 _ 


KGRAO 


KPLOT 


CODE  FOR  REFRACTIVITY 

_ PROFILE  PRINTOUT. _ _ _ 

<PRO=0  NO  PRINTOUT 

_  <PR0=1  PRINTOUT  _ 

CODE  FOR  REFRACTIVITY  GRADIENT 

_ PROFILE  PRINTOUT. _ 

<GRAD=0  NO  PRINTOUT 

<GRAD=1  PRINTOUT  _ 

**  HIS  VARIABLE  IS  NOT  USED. 

_  _ SET  KRAY=Q _ _ _ 

CODE  FOR  RELATIVE  FIELO  STRENGTH 

_OR  POWER  DENSITY  PRINTOUT. _ 

KPLOT=0  MO  PRINTOUT 
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TABLE  10  (continued) 


<PL0T=1 


<PLOT  =  2 


RELATIVE  FIE .P 
STRENGTH  PRINTOUT 
RELATIVE  POWER 
DENSITY  PRINTOUT  . 


NOATA 
♦  10  _ 


NGRAD 


CODE  FOR  REFRACTIVITY 

_ PROFILE  PLOT. _ 

NPRO=Q  NO  PLOT 

_ _ »PR0=1  PLOT _ _ 

CODE  FOR  REFRACTIVITY  GRADIENT 

_ PROFILE  PLOT. _ 

MGRAD=Q  NO  PLOT 

_ NGRAD=1 .  PLOT  _ 

CODE  FOR  RAY  TRACE  PLOT. 

_  NRAY=Q  NO  PLOT _ 

NRAY=1  PLOT 

_ CODE  FOR  RELATIVE  FIELO. STRENGTH _ 

OR  POWER  DENSITY  PLOT. 

_ NPLOT=D..  NO  PLOT _ 

NPL0T=1  RELATIVE  FIELO 

_ _  STRENGTH  PLOT _ 

N PLOT* 2  RELATIVE  POWER 
_ DEN  S I  T_Y_PL  01 _ 


NPL0T*2 


NOATA  SCALE 

..+11. _ 


F15.7  SCALE  FACTOR  FOR  ENLARGING  OR 

_ REDUCING  PLOT  SIZE  FROM  ITS _ 

NORMAL  5  X  10  INCH  FORMAT. 

_ NORMAL  PLOT  SIZE  GIVEN  HITH 

SCALE*1.G 

_ ♦*  (OMIT  THIS  CARO  AND  ALL  CAROS 

THAT  FOLLOW  IF  THERE  ARE  NO 
_ PLO  T  SI _ 


NOATA 

IAHO 

12 

NUMBER  OF  CHARACTERS  IN  PLOT 

♦  12 _ 

BANNER. 

AHD 

7A13 

CHARACTERS  IN  PLOT  BANNER. 

••  (THIS  CARD  IS  ALWAYS  REAO 

WHEN  THERE  ARE  ANY  PLOTS* 

IX.TL _ 


NUMBER  OF  CHARACTERS  IN  X-AXIS__ 
LABEL. 

CHARACTERS  IN  X-AXIS  LABEL.  _ 

••(THIS  CARO  READ  FOR  ALL  PLOTS* 


NOATA 
♦  1<* 


NUMBER  OF  CHARACTERS  IN.  Y-AXIS _ 

LABEL. 

CHARACTERS  IN  Y-AXIS  LABEL. 
'••(THIS  CARD  READ  FOR  ALL  PLOTS* 


NOATA 
♦  15 


number  off  Characters  in  first 

LINE  OF  PLOT  TITLE. 
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TABLE  10  (continued) 


TTL  7A10  CHARACTERS  IN  FIRST  LINE  OF  PLOT 

_  TITLE.  _  __ 

** (THIS  CARO  READ  FOR  ALL  PLOTS) 


NOATA _ ILLS _ _  12 _ NUMBER.  OF  CHARACTERS  _IN_SECONO 

♦16  LINE  OF  PLOT  TITLE. 

_ T.LI _ 7413 _ CHARACTERS  IN  SECOND  LINE  0F_ 

PLOT  TITLE. 

_ **_miS  CARO  REAO  FOR  ALL  PLOTS) 


*****  REPEAT  CAROS  (NOATA+ll)  THROUGH  (NOATA+16)  FOR  EACH  PLOT 
_ IN  _A  GIVEN  SIMULATION  .RUN. _ 

*****  REPEAT  CAR0S_(2 )_T_H ROUGH  (N0ATA4-16)  _FOR  EACH  SIMULATION 
RUN. 
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TABLE  12  (Continued) 
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